{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import peregrine\n",
    "import sys\n",
    "from peregrine._falcon4py import ffi as falcon_ffi\n",
    "from peregrine._falcon4py import lib as falcon4py\n",
    "from peregrine._shimmer4py import ffi as shimmer_ffi\n",
    "from peregrine._shimmer4py import lib as shimmer4py\n",
    "import numpy as np\n",
    "import mmap\n",
    "from pprint import pprint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "f=open(\"/mnt/data/test3_28x_wd_2/wd/index/seq_dataset.seqdb\", \"rb\")\n",
    "seqdb=mmap.mmap(f.fileno(), 0, flags=mmap.MAP_SHARED, prot=mmap.PROT_READ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmap = dict(zip(b\"ACGT\", b\"TGCA\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "read_data = {}\n",
    "with open(\"/mnt/data/test3_28x_wd_2/wd/index/seq_dataset.idx\") as f:\n",
    "    for row in f:\n",
    "        row = row.strip().split()\n",
    "        rid, rname, rlen, offset = row\n",
    "        rid = int(rid)\n",
    "        rlen = int(rlen)\n",
    "        offset = int(offset)\n",
    "        read_data.setdefault(rid, {})\n",
    "        read_data[rid][\"name\"] = rname\n",
    "        read_data[rid][\"length\"] = rlen\n",
    "        read_data[rid][\"offset\"] = offset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "py_mmer_L2=shimmer_ffi.new(\"py_mmer_t *\")\n",
    "shimmer4py.build_shimmer_map4py(py_mmer_L2, \n",
    "                             b\"/mnt/data/test3_28x_wd_2/wd/index/seq_dataset\", \n",
    "                             b\"/mnt/data/test3_28x_wd_2/wd/index/shimmer-L2\", \n",
    "                             1, 1, 2, 240)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the range of the shimmers of each in the py_mmer vector\n",
    "def get_mmer_index(py_mmer):\n",
    "    mmer_index = {}\n",
    "    for i in range(py_mmer.mmers.n):\n",
    "        mmer = py_mmer.mmers.a[i]\n",
    "        rid = mmer.y >> 32\n",
    "        mmer_index.setdefault(rid, [None, None])\n",
    "        if mmer_index[rid][0] == None:\n",
    "            mmer_index[rid][0] = i\n",
    "        if mmer_index[rid][1] == None or i > mmer_index[rid][1]:\n",
    "            mmer_index[rid][1] = i\n",
    "    return mmer_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "mmer_index_L2 = get_mmer_index(py_mmer_L2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the range of shimmers for the read 3 is from 18261594 to 18261629\n"
     ]
    }
   ],
   "source": [
    "rid = 3\n",
    "s,e=mmer_index_L2[rid]\n",
    "print(f\"the range of shimmers for the read {rid} is from {s} to {e}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "basemap = {1:\"A\",2:\"C\",4:\"G\",8:\"T\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_shimmers_for_read(rid, mmer_index, py_mmer, mc_l=2, mc_h=240):\n",
    "    kset = set()\n",
    "    shimmers = []\n",
    "    s,e = mmer_index[rid]\n",
    "    for i in range(s,e):\n",
    "        x = py_mmer.mmers.a[i].x\n",
    "        y = py_mmer.mmers.a[i].y\n",
    "        span = x & 0xFF\n",
    "        mmer = x >> 8\n",
    "        mmcount = shimmer4py.get_mmer_count(py_mmer, mmer)\n",
    "        if mmcount < mc_l or mmcount > mc_h:\n",
    "            continue\n",
    "        rid = y >> 32\n",
    "        pos_end = ((y & 0xFFFFFFFF) >> 1) + 1\n",
    "        strand = y & 0x1\n",
    "        \n",
    "        mm_str = \"{:014X}\".format(mmer)\n",
    "        \n",
    "        s = read_data[rid][\"offset\"]\n",
    "        e = s + read_data[rid][\"length\"]\n",
    "        kmer =  \"\".join([basemap[c&0x0F] for c in seqdb[s:e][pos_end-span:pos_end]])\n",
    "        \n",
    "        shimmers.append( (rid, span, pos_end, strand, mmer, mm_str, kmer, mmcount) )\n",
    "        kset.add(mmer)\n",
    "    return shimmers, kset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "rid = 1\n",
    "shimmers_L2, kset_L2 = get_shimmers_for_read(rid, mmer_index_L2, py_mmer_L2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L2 shimmer for read 1\n",
      "19\n"
     ]
    }
   ],
   "source": [
    "print(f\"L2 shimmer for read {rid}\")\n",
    "pprint(len(shimmers_L2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "shimmer_count_L2 = []\n",
    "all_rids = list(read_data.keys())\n",
    "for rid in all_rids[:50000]:\n",
    "    if rid not in mmer_index_L2:\n",
    "        continue\n",
    "    shimmers_L2, kset_L2 = get_shimmers_for_read(rid, mmer_index_L2, py_mmer_L2)\n",
    "    if len(shimmers_L2) < 2:\n",
    "           continue\n",
    "    shimmer_count_L2.append( (len(shimmers_L2),read_data[rid][\"length\"]/(len(shimmers_L2)-1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean distance: 819.3684715811905\n",
      "median distance: 633.6\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcHFW5//HPl7CTQEIIGpJAEFABRZaw6VXx6mUTQQU0XJSA0agXVER+LBcVVFSUiyDXNUDY9BoQBaKigAiCSICEHYImhkAmCSTsiSwm8Pz+OKeZyqR7pisz3T3L9/16zWuqT52uemrpeqrq1KKIwMzMrF5rtDoAMzPrW5w4zMysFCcOMzMrxYnDzMxKceIwM7NSnDjMzKwUJ44mkPQTSV/poWFtLmmZpEH5802SPtkTw87D+72kCT01vBLjPV3Sk5Ier7P+aZJ+1qBYjpT0l8LnZZLe0IhxNZqkiySd3uo4mq3jMqzSvyXreX/hxNFNkuZJelHSUknPSvqrpM9Iem3eRsRnIuIbdQ7rfZ3ViYjHImJwRLzSA7GvsvGNiP0i4uLuDrtkHGOALwHbRcTrq/TfS1JbM2MqyvN7bmd1Wh1jI/TnpNOK9bw/ceLoGR+IiCHAFsAZwInABT09Eklr9vQwe4ktgKciYnGrAzFrpT7zG48I/3XjD5gHvK9D2W7Aq8Bb8ueLgNNz9ybAb4FngaeBW0gJ/NL8nReBZcAJwFgggInAY8DNhbI18/BuAr4N3AE8B1wNbJz77QW0VYsX2Bf4F7A8j+/ewvA+mbvXAL4MPAosBi4BNsr9KnFMyLE9CZzSyXzaKH9/SR7el/Pw35en+dUcx0UdvrdBh/7LgM2A04DL8zCXAg8C4wrf2wz4VR7fI8DnO4ltODANeD7Px28Afyn0D2Dr3L0/8FAe5wLg+E5i3A24LS/rRcAPgLU7DPczwGzgGeCHgAr9PwXMyuN6CNh5NabtIuAnwPV5OH8Gtij0f3Pu9zTwN+AjuXxSXjf+lafnN8BRwG8K350DXF74PB/YsbPh5n7rAP+T15sncnzrFddZ0hHo4jzfjupk+o4E5uZpewQ4vFD+lzyeZ3K//Qrfu4n29fxI4Fbg7Lys5gJvz+XzcxwTOszTHwG/z/PmVuD1wDl5XA8DO9WzLpLW4yuAn5HWv0+S1psZ+fMTwPdavZ1bZb63OoC+/keVxJHLHwM+W1jRKonj2/mHslb+eyd5Y9FxWLRvnC8hbZzWo3riWAC8Jdf5FfCz3G8vaiSO3H1apW6hf/EH9QnSxuENwGDg18ClHWI7L8f1NuBlYNsa8+kSUlIbkr/7d2BirTg7fLfadJwGvETakA/K83V67rcGMBP4KrB2jn8usE+N4U8lJaEN8nxcQO3EsQh4Z+4eRvvGvFqMuwB7AGvmaZ4FHNthuL8FhgKbkzYs++Z+h+Y4dgUEbE06Mis7bReRNqrvIm2wv1+Ztjy980kJYU1gZ9IOwPYd19v8+Q2kDesawEjSDsCCQr9ncr+uhnsOKVFvnNeH3wDfLszHFcDXSb+P/YEXgGFVpm0D0sb1TfnzyMI4jiQlvk/l9eOzwELaf2s3sXLiWJHjHQScTvr9/jDPs73zPBxcmC9P5uW7LvAnUkI4ovD9G+tZF0nr8XLgg7nueqSdjY/n/oOBPVq9nev451NVjbOQ9MPoaDlpBd8iIpZHxC2R15BOnBYR/4yIF2v0vzQiHoiIfwJfAT5SaTzvpsNJeztzI2IZcDIwvsPh9Nci4sWIuBe4l5RAVpJj+ShwckQsjYh5wFnAx7sZ318i4ppI7T2XFsa9KzAiIr4eEf+K1D5xHjC+RmwHA1/N8/gBoLNz38uB7SRtGBHPRMRdtSpGxMyImB4RK/I0/xR4d4dqZ0TEsxHxGHAjsGMu/yTw3Yi4M5I5EfFomWkr+F1E3BwRLwOnAHvmdqUDgHkRcWGO8S7SjschNaansme/Y56Oa4EFkt6cP98SEa92NlxJIm3MvxgRT0fEUuBbHeJfDnw9/z6uIe3Vv6nGtL0KvEXSehGxKCIeLPR7NCLOy+vHxaTf3etqDOeRHO8rwGXAmBzDyxFxHenIa+tC/Svz8n0JuBJ4KSIuKXx/p1yvnuV1W0RcFRGv5t/4cmBrSZtExLKImF4j5pZx4micUaTD9I7OJO3FXydprqST6hjW/BL9HyXtqW1SV5Sd2ywPrzjsNVn5x1e8CuoF0h5SR5uQ9rY6DmtUN+PrOO51c1LbAtgsX6zwrKRngf+m+kZjBGmaOs7DWg4m7QU/KunPkvasVVHSGyX9VtLjkp4nbSA7Lpda828M8I8qgy0zbRWvTVveAXiatGy3AHbvMKzDSaddavkz6ajgXbn7JlLSeHf+XImx1nBHAOsDMwv9/pDLK56KiBWFz1XXq7yj9FHS6b5Fkn6Xk1jF44W6L+TOausnpFNCFS/m73QsG9xJ/Vp161leHX/fE4E3Ag9LulPSATVibpm+0RDTx0jalbRRXOVywLyH9SXgS5K2B26UdGdE3EA6dVFNV0ckYwrdm5P2WJ4E/kn6kVbiGsTKP9CuhruQtOIXh72C9CMZ3cV3i57MMW1BOldfGdaCOr9f9hHO80l7kNvUUXcJaZrGkM5NV2KrHkjEncBBktYCjiGd4hpTI8YfA3cDh0XEUknHUmNvvsY0bFWjvN5pq3ht/ZA0mHQkvDAP688R8R81vldtmv4MfADYkpQIK0lhT1IbTiXGqsPNVxu+SDqlVO/yrykirgWulbQe6RTReaTTv71FPctrpfkcEbOBw/K8+jBwhaThOVH2Cj7i6EGSNsx7B1NJbQf3V6lzgKSt8yH788Ar+Q/SBnl17hf4mKTtJK1POjd8RT5k/jtpL/z9eUP3ZdI524ongLHFS4c7+AXwRUlb5g3Ot4DLOuwNdinHcjnwTUlDJG0BHEdqEKzHE8BwSRvVWf8O4HlJJ0paT9IgSW/JCb1abL8GTpO0vqTtSA3+q5C0tqTDJW0UEctpX361YhyS6yzLe8KfrTN+gPOB4yXtomTrPN/qnraC/SX9m6S1SQ3/t0fEfFL7yhslfVzSWvlvV0nbFqap4/r4Z+A9pMbsNtLFHfuSLjC4O9epOdx8Kus84GxJm+b5OkrSPiXmDfl7r5N0oKQNSO1ry2hfHr1F6eUl6WOSRuR59Wwu7lXT5cTRM34jaSlp7+IU4HukhrZqtgH+SFrJbwN+FBE35X7fBr6cD2mPLzH+S0kNdo+TGus+DxARzwH/RdoILSAdgRTvNfhl/v+UpGrn6qfkYd9Mavx7CfhcibiKPpfHP5d0JPZ/efhdioiHSUlsbp43m3VR/xXSXvGOOe4nSfOgVuI5hnRq4XHSfLywk8F/HJiXTz19BvhYJzEeD/wnqV3gPNK577pExC+Bb5Lm01LgKtLVcmWnjTyMU0mnqHYhHSFUjn73Jp1vX5in/zu071xcQGrPeVbSVfk7fyetu7fkz8+TlumtObZ6hnsi6XTt9Dwf/0jtNozOrEE6el+Yp+3dpPW911jN5bUv8KCkZaSLGcbntpReo3KFgZmZWV18xGFmZqU4cZiZWSlOHGZmVooTh5mZldIv7+PYZJNNYuzYsa0Ow8ysT5k5c+aTETGiq3r9MnGMHTuWGTNmtDoMM7M+RVJnT014jU9VmZlZKU4cZmZWihOHmZmV0i/bOMzMipYvX05bWxsvvdSrntzRMuuuuy6jR49mrbXWWq3vO3GYWb/X1tbGkCFDGDt2LOn5ogNXRPDUU0/R1tbGlltuuVrD8KkqM+v3XnrpJYYPHz7gkwaAJIYPH96toy8nDjMbEJw02nV3XjhxmJlZKW7jMLMBZ/zk23p0eFMn1XyD8GsGDx7MsmXLenS81YY5f/58jjjiCB5//HHWWGMNJk2axBe+8IUeHW/DjjgkTZG0WNIDVfodLykkbZI/S9K5kuZIuk/SzoW6EyTNzn9V38xmvV9P/1DNrLo111yTs846i1mzZjF9+nR++MMf8tBDD3X9xRIaearqItKbrFYiaQzwH8BjheL9SG/G2waYRHpXM5I2Jr25bHdgN+BUScMaGLOZWcOdeeaZ7Lrrruywww6ceuqpAJx44on86Ec/eq3OaaedxllnnVWzfi0jR45k553TvveQIUPYdtttWbCg2693X0nDEkdE3Ex6nWNHZwMnsPIL2g8CLolkOjBU0khgH+D6iHg6Ip4BrqdKMjIz6yuuu+46Zs+ezR133ME999zDzJkzufnmmxk/fjyXXdb+duHLL7+cQw89tGb9esybN4+7776b3XffvUenoaltHJIOBBZExL0dWvVHkd7XXdGWy2qVVxv2JNLRCptvvnkPRm1m1nOuu+46rrvuOnbaaScAli1bxuzZs5k4cSKLFy9m4cKFLFmyhGHDhrH55ptz7rnnVq3/rne9q9PxLFu2jIMPPphzzjmHDTfcsEenoWmJQ9L6wCmkl9iv0rtKWXRSvmphxGRgMsC4ceP8InUz65UigpNPPplPf/rTq/Q75JBDuOKKK3j88ccZP358l/VrWb58OQcffDCHH344H/7wh3ss9opmXo67FbAlcK+kecBo4C5JrycdSYwp1B0NLOyk3MysT9pnn32YMmXKa1dDLViwgMWLFwMwfvx4pk6dyhVXXMEhhxzSZf1qIoKJEyey7bbbctxxxzVkGpp2xBER9wObVj7n5DEuIp6UNA04RtJUUkP4cxGxSNK1wLcKDeJ7Ayc3K2Yz65/quXy2Ufbee29mzZrFnnumGAYPHszPfvYzNt10U7bffnuWLl3KqFGjGDlyZJf1q7n11lu59NJLeetb38qOO+4IwLe+9S3233//HpsGRTTmrI6kXwB7AZsATwCnRsQFhf7zaE8cAn5Aavh+ATgqImbkep8A/jt/7ZsRcWFX4x43blz4RU69y/jJt7X0x2oD26xZs9h2221bHUavUm2eSJoZEeO6+m7Djjgi4rAu+o8tdAdwdI16U4ApPRqcmZmtNj9yxMzMSnHiMLMBoVGn5fui7s4LJw5rmvGTb/OjR6wl1l13XZ566iknD9rfx7Huuuuu9jD8kEMz6/dGjx5NW1sbS5YsaXUovULlDYCry4nDzPq9tdZaa7Xfdmer8qkqMzMrxYnDms7tHGZ9mxOHmZmV4jYOaygfXZj1Pz7iMDOzUpw4zMysFCcOaxifpjLrn5w4rCV8F7lZ3+XEYWZmpThxmJlZKU4c1qN8+sms//N9HNYjnDDMBg4fcZiZWSlOHGZmVopPVVlLFU9xTZ20ZwsjMbN6NeyIQ9IUSYslPVAoO1PSw5Luk3SlpKGFfidLmiPpb5L2KZTvm8vmSDqpUfFaz/E9Gmb9WyNPVV0E7Nuh7HrgLRGxA/B34GQASdsB44Ht83d+JGmQpEHAD4H9gO2Aw3JdMzNrkYYljoi4GXi6Q9l1EbEif5wOVN5deBAwNSJejohHgDnAbvlvTkTMjYh/AVNzXTMza5FWNo5/Avh97h4FzC/0a8tltcpXIWmSpBmSZvi9wn2bT3OZ9W4tSRySTgFWAD+vFFWpFp2Ur1oYMTkixkXEuBEjRvRMoGZmtoqmX1UlaQJwAPDeiKgkgTZgTKHaaGBh7q5VbmZmLdDUIw5J+wInAgdGxAuFXtOA8ZLWkbQlsA1wB3AnsI2kLSWtTWpAn9bMmM3MbGUNO+KQ9AtgL2ATSW3AqaSrqNYBrpcEMD0iPhMRD0q6HHiIdArr6Ih4JQ/nGOBaYBAwJSIebFTMZmbWtYYljog4rErxBZ3U/ybwzSrl1wDX9GBo1ku5Udysb/AjR8zMrBQnDjMzK8WJw8zMSnHiMDOzUpw4rNvcqG02sDhxmJlZKU4cZmZWihOHmZmV4sRhZmalOHGYmVkpThzWq/mKLbPep+mPVTerhxOGWe/lIw4zMyvFRxy22nxUYDYw+YjDzMxKceIwM7NSnDjMzKwUJw4zMyvFicPMzEpx4jAzs1IaljgkTZG0WNIDhbKNJV0vaXb+PyyXS9K5kuZIuk/SzoXvTMj1Z0ua0Kh4rRxfims2cDXyiOMiYN8OZScBN0TENsAN+TPAfsA2+W8S8GNIiQY4Fdgd2A04tZJsbOBxsjLrHRqWOCLiZuDpDsUHARfn7ouBDxbKL4lkOjBU0khgH+D6iHg6Ip4BrmfVZGRmZk3U7DvHXxcRiwAiYpGkTXP5KGB+oV5bLqtVvgpJk0hHK2y++eY9HLa1ko80zHqX3tI4ripl0Un5qoURkyNiXESMGzFiRI8GZ2Zm7ZqdOJ7Ip6DI/xfn8jZgTKHeaGBhJ+VmZtYizU4c04DKlVETgKsL5Ufkq6v2AJ7Lp7SuBfaWNCw3iu+dy8zMrEUa1sYh6RfAXsAmktpIV0edAVwuaSLwGHBorn4NsD8wB3gBOAogIp6W9A3gzlzv6xHRscHdzMyaqGGJIyIOq9HrvVXqBnB0jeFMAab0YGhmZtYNfh+HleIrnMyst1xVZWZmfYQTh5mZleLEYWZmpThxmJlZKU4cZmZWSpeJQ9Kl9ZSZNZOv7jJrnXqOOLYvfpA0CNilMeGYmVlvVzNxSDpZ0lJgB0nP57+lpOdLXV3re2Zm1r/VTBwR8e2IGAKcGREb5r8hETE8Ik5uYoxmZtaLdHnneEScLGkUsEWxfn5Rk1lTuW3DrPW6TBySzgDGAw8Br+TiAJw4zMwGoHqeVfUh4E0R8XKjgzEzs96vnquq5gJrNToQMzPrG+o54ngBuEfSDcBrRx0R8fmGRWVmZr1WPYljWv4zMzOr66qqi5sRiJmZ9Q31XFX1COkqqpVExBsaEpGZmfVq9ZyqGlfoXpf0nvCNGxOOmZn1dl1eVRURTxX+FkTEOcC/NyE2MzPrheo5VbVz4eMapCOQId0ZqaQvAp8knQK7HzgKGAlMJR3N3AV8PCL+JWkd4BLSgxWfAj4aEfO6M34zM1t99ZyqOqvQvQKYB3xkdUeYH1/yeWC7iHhR0uWkO9P3B86OiKmSfgJMBH6c/z8TEVtLGg98B/jo6o7fzMy6p56rqt7ToPGuJ2k5sD6wiHT66z9z/4uB00iJ46DcDXAF8ANJiohVGuzNzKzx6nmR00aSvidpRv47S9JGqzvCiFgA/A/wGClhPAfMBJ6NiBW5WhswKnePAubn767I9YdXiXNSJcYlS5asbnjWh4yffJsfemjWAvU8cmQKsJR0euojwPPAhas7QknDSEcRWwKbARsA+1WpWjmiUCf92gsiJkfEuIgYN2LEiNUNz8zMulBP4tgqIk6NiLn572tAd+7heB/wSEQsiYjlwK+BtwNDJVVOnY0GFubuNmAMQO6/EfB0N8Zv/YyPOsyaq57E8aKkf6t8kPQO4MVujPMxYA9J60sS8F7SI9tvBA7JdSbQ/pbBafkzuf+f3L5hZtY69VxV9Vng4kK7xjPAkas7woi4XdIVpEtuVwB3A5OB3wFTJZ2eyy7IX7kAuFTSHNKRxvjVHbeZmXVfPVdV3QO8TdKG+fPz3R1pRJwKnNqheC6wW5W6L5HuVrcW8ykhM4P6rqr6lqShEfF8RDwvaVg+KjAzswGonjaO/SLi2cqHiHiGdLOemZkNQPUkjkH5sR8ASFoPWKeT+mZm1o/V0zj+M+AGSReS7p/4BOnObrNeo9L+MnXSni2OxKz/q6dx/LuS7iPdfyHgGxFxbcMjMzOzXqmeIw4i4g/AHxoci5mZ9QH1tHGY9Sl+hpVZYzlxmJlZKU4cZmZWSs02Dkn3U+UptBURsUNDIjIzs16ts8bxA/L/o/P/S/P/w4EXGhaR9TpuLzCzopqJIyIehfQ03Ih4R6HXSZJuBb7e6ODMynKSM2u8eto4NujwWPW3k16+ZGZmA1A993FMBKYUHqv+LOnucTMzG4DquXN8Ju2PVVdEPNf4sMzMrLeq685xSe8HtgfWTS/tg4hwG4eZ2QBUz/s4fgJ8FPgc6VlVhwJbNDguMzPrpeppHH97RBwBPBMRXwP2BMY0NiwzM+ut6kkcL+b/L0jaDFgObNm4kMzMrDerp43jt5KGAmcCd5HuJj+/oVGZmVmv1eURR0R8IyKejYhfkdo23hwRX+nOSCUNlXSFpIclzZK0p6SNJV0vaXb+PyzXlaRzJc2RdJ+knbszbhs4fDOgWWPU0zi+vqSvSDovIl4GNpV0QFff68L3gT9ExJuBtwGzgJOAGyJiG+CG/BlgP2Cb/DcJ+HE3x21mZt1QTxvHhcDLpEZxgDbg9NUdYb4f5F3ABQAR8a+IeBY4iPZX0l4MfDB3HwRcEsl0YKikkas7fivHe+1m1lE9iWOriPguqVGciHiRdFnu6noDsAS4UNLdks6XtAHwuohYlMexCNg01x8FzC98vy2XrUTSJEkzJM1YsmRJN8IzM7PO1JM4/iVpPfIj1iVtRToCWV1rAjsDP46InYB/0n5aqppqSWqVx71HxOSIGBcR40aMGNGN8MzMrDP1JI5TSe8bHyPp56T2hxO6Mc42oC0ibs+fryAlkicqp6Dy/8WF+sX7RkYDC7sxfjMz64ZOE4fS80UeBj4MHAn8AhgXETet7ggj4nFgvqQ35aL3Ag8B04AJuWwCcHXungYcka+u2gN4rnJKy8zMmq/T+zgiIiRdFRG7AL/rwfF+Dvi5pLWBucBRpCR2uaSJwGOkR5sAXAPsD8whvUDqqB6Mw8zMSqrnBsDpknaNiDt7aqQRcQ8wrkqv91apG7S/hdCslMpVYVMn7dlFTTOrVz1tHO8BbpP0j3wD3v2S7mt0YGY9yZcVm/Wceo449mt4FGZN4KMPs55Rz4ucHm1GIGZm1jfUc6rKzMzsNU4cZmZWihOHmZmV4sRhZmalOHGYmVkpThxmZlaKE4cNOL4Z0Kx76rkB0AYgb1zNrBYfcZiZWSlOHGZmVooTh5mZleLEYWZmpbhx3AakYuO/n5ZrVo6POMzMrBQnDjMzK8WJw1bhezjMrDNOHGZmVkrLEoekQZLulvTb/HlLSbdLmi3pMklr5/J18uc5uf/YVsVsZmatPeL4AjCr8Pk7wNkRsQ3wDDAxl08EnomIrYGzcz0zM2uRliQOSaOB9wPn588C/h24Ile5GPhg7j4ofyb3f2+ub2ZmLdCqI45zgBOAV/Pn4cCzEbEif24DRuXuUcB8gNz/uVzfrEeMn3ybLwgwK6HpiUPSAcDiiJhZLK5SNeroVxzuJEkzJM1YsmRJD0RqZmbVtOKI4x3AgZLmAVNJp6jOAYZKqtzJPhpYmLvbgDEAuf9GwNMdBxoRkyNiXESMGzFiRGOnwPo1H32Yda7piSMiTo6I0RExFhgP/CkiDgduBA7J1SYAV+fuafkzuf+fImKVIw4zM2uO3nQfx4nAcZLmkNowLsjlFwDDc/lxwEktis/MzGjxQw4j4ibgptw9F9itSp2XgEObGpiZmdXkp+OaZW7bMKuPE4dZFX7sulltvamNw8zM+gAfcdhrfKrGzOrhIw4zMyvFicPMzEpx4jAzs1KcOMzMrBQnDrMu+Om5Zitz4jAzs1KcOMzMrBQnDjMzK8WJw2w1uM3DBjInDjMzK8WPHDErwUcaZj7iMKubk4ZZ4sRhtpp8f4cNVE4cBnhv2szq58RhZmalOHGYdZOP1mygceIw60FOIjYQND1xSBoj6UZJsyQ9KOkLuXxjSddLmp3/D8vlknSupDmS7pO0c7NjNjOzdq044lgBfCkitgX2AI6WtB1wEnBDRGwD3JA/A+wHbJP/JgE/bn7IZmZW0fQbACNiEbAody+VNAsYBRwE7JWrXQzcBJyYyy+JiACmSxoqaWQejlmvUOsUVaV86qQ9mxmOWUO19M5xSWOBnYDbgddVkkFELJK0aa42Cphf+FpbLlspcUiaRDoiYfPNN29o3GadcTuH9XctSxySBgO/Ao6NiOcl1axapSxWKYiYDEwGGDdu3Cr9rTpv5MysrJZcVSVpLVLS+HlE/DoXPyFpZO4/Elicy9uAMYWvjwYWNitWs57iO82tv2jFVVUCLgBmRcT3Cr2mARNy9wTg6kL5Efnqqj2A59y+YWbWOq04VfUO4OPA/ZLuyWX/DZwBXC5pIvAYcGjudw2wPzAHeAE4qrnhmnWfjzSsP2nFVVV/oXq7BcB7q9QP4OiGBmVmZnXzneNmZlaKE4dZk/m0lfV1ThxmLeArrKwvc+IYoLzh6h28DKwvcuIwM7NSWvrIETNrVzz68LOtrDdz4jBrsWqnq4oPR3RCsd7Gp6rMzKwUJw6zXqzj0YgvarDewInDzMxKceIYgLzHOjB4OVujOHGY9WM+tWWN4MRh1se5HcSaTenhs/3LuHHjYsaMGa0Oo9fyRmVg8yW9VoukmRExrqt6PuIws7p4h8MqnDjMBpjuJIBqp8GcUAYe3zk+gPgHbhVl7kyvdWe7T3kNXG7jGCCcNKwZionIiaXvcRuHmTVd2Su8Kv19+qtvceIws4Yobvw7dtdKKPWWNzqxOHF1rs+cqpK0L/B9YBBwfkScUauuT1W18w/A+pPK6a961uuOp8rqOYVWbbi1Tr/VKuvLp+jqPVXVJxKHpEHA34H/ANqAO4HDIuKhavWdOBInDRvoyiSaZo2/bGIpJqNaiamnHr3f3xLHnsBpEbFP/nwyQER8u1r9gZw4nCzMBrZmJI6+cjnuKGB+4XMbsHuxgqRJwKT8cZmkv3VjfJsAT3bj+43iuMpxXOU4rnJ6ZVyXfbpbcW1RT6W+kjhUpWylQ6WImAxM7pGRSTPqybrN5rjKcVzlOK5yBnJcfeWqqjZgTOHzaGBhi2IxMxvQ+kriuBPYRtKWktYGxgPTWhyTmdmA1CdOVUXECknHANeSLsedEhEPNnCUPXLKqwEcVzmOqxzHVc6AjatPXFVlZma9R185VWVmZr2EE4eZmZXixFEgaV9Jf5M0R9JJTR73GEk3Spol6UFJX8jlp0laIOme/Ld/4Tsn51j/JmmfBsY2T9L9efwzctnGkq6XNDv/H5bLJencHNd9knZuUExvKsyTeyQ9L+nYVswvSVMkLZb0QKGs9PyRNCHXny1pQoPiOlPSw3ncV0oamsvHSnqxMN9+UvjOLnn5z8mxV7s8vrtxlV5uPf17rRHXZYWY5km6J5c3c37V2ja0bh2LCP+ldp4nNO8OAAALhElEQVRBwD+ANwBrA/cC2zVx/COBnXP3ENIjVrYDTgOOr1J/uxzjOsCWOfZBDYptHrBJh7LvAifl7pOA7+Tu/YHfk+692QO4vUnL7nHSzUtNn1/Au4CdgQdWd/4AGwNz8/9huXtYA+LaG1gzd3+nENfYYr0Ow7kD2DPH/HtgvwbEVWq5NeL3Wi2uDv3PAr7agvlVa9vQsnXMRxztdgPmRMTciPgXMBU4qFkjj4hFEXFX7l4KzCLdMV/LQcDUiHg5Ih4B5pCmoVkOAi7O3RcDHyyUXxLJdGCopJENjuW9wD8i4tFO6jRsfkXEzcDTVcZXZv7sA1wfEU9HxDPA9cC+PR1XRFwXESvyx+mke6JqyrFtGBG3Rdr6XFKYlh6LqxO1lluP/147iysfNXwE+EVnw2jQ/Kq1bWjZOubE0a7aY00623A3jKSxwE7A7bnomHzIOaVyOEpz4w3gOkkzlR7tAvC6iFgEacUGNm1BXBXjWfkH3er5BeXnTyvm2ydIe6YVW0q6W9KfJb0zl43KsTQjrjLLrdnz653AExExu1DW9PnVYdvQsnXMiaNdl481aUoQ0mDgV8CxEfE88GNgK2BHYBHpcBmaG+87ImJnYD/gaEnv6qRuU+ej0g2hBwK/zEW9YX51plYczZ5vpwArgJ/nokXA5hGxE3Ac8H+SNmxiXGWXW7OX52GsvHPS9PlVZdtQs2qNGHosNieOdi1/rImktUgrxs8j4tcAEfFERLwSEa8C59F+eqVp8UbEwvx/MXBljuGJyimo/H9xs+PK9gPuiogncowtn19Z2fnTtPhyo+gBwOH5dAr5VNBTuXsmqf3gjTmu4umshsS1GsutmfNrTeDDwGWFeJs6v6ptG2jhOubE0a6ljzXJ51AvAGZFxPcK5cX2gQ8BlSs+pgHjJa0jaUtgG1KjXE/HtYGkIZVuUuPqA3n8lasyJgBXF+I6Il/ZsQfwXOVwukFW2hNs9fwqKDt/rgX2ljQsn6bZO5f1KKUXop0IHBgRLxTKRyi99wZJbyDNn7k5tqWS9sjr6BGFaenJuMout2b+Xt8HPBwRr52Caub8qrVtoJXrWHda+/vbH+lqhL+T9h5OafK4/4102HgfcE/+2x+4FLg/l08DRha+c0qO9W9088qNTuJ6A+mKlXuBByvzBRgO3ADMzv83zuUCfpjjuh8Y18B5tj7wFLBRoazp84uUuBYBy0l7dRNXZ/6Q2hzm5L+jGhTXHNJ57so69pNc9+C8fO8F7gI+UBjOONKG/B/AD8hPnOjhuEovt57+vVaLK5dfBHymQ91mzq9a24aWrWN+5IiZmZXiU1VmZlaKE4eZmZXixGFmZqU4cZiZWSlOHGZmVooTh/UKSk9HPT53f13S+zqp+0FJ2zUvupXG/VqcddYfKum/GhlTT+pq3uc6B6rk02glXaP8JN5Gk7STpPNzd6nllb/zx8IjT6wKJw7rVOUmp2aKiK9GxB87qfJB0tNB+4KhQJ9JHHXMeyJiWkScUXK4+0fEs92JrcS6+N/A/3ZjVJfSh5ZZKzhx9COSrsoPInyw8jBCSZ+V9N1CnSMl/W/u/pikO5TeJ/DTwp2wy/Ke5+3AnpK+KulOSQ9ImpzvZEXSrvmhdLcpvefhgVw+KH++M/f/dI14T1F6n8IfgTcVyi+SdEjuPkPSQ3k4/yPp7aRnU52Z495K0qfyuO6V9CtJ6xeGc66kv0qaWxlm7neC0jsT7pV0Ri7bStIf8jy8RdKba8zqt0n6k9I7DT5VGOb/K0zz13LxGcBWOdYzJf1I0oG5/pWSpuTuiZJO72K57J3n9V2Sfqn07KLK+1K+lsvvrxZ3Xu5XSfqNpEckHSPpOKWH9E2XtHGVeV91uHlYPyjU/7HS+yLmSnq30kMKZ0m6qDD+eZI2kfQZtb/D4hFJN9YxbV+V9BfgUEmfL6wPU6tM5xBgh4i4t7PlJWkvSTfnZfCQpJ9IqmwPp5GeSGC19NTds/5r/R/td46uR7pzdTgwgvT46Uqd35PuRN0W+A2wVi7/EXBE7g7gIx2Hm7svJd8lm8fx9tx9Bvn9BMAk4Mu5ex1gBrBlh1h3Id3Vuj6wIelO1uNzv4uAQ0jvDfgbvHaj6tBi/8Kwhhe6Twc+V6j3S9IO0naV+UB6vtVfgfU7zLcbgG1y9+7An6rM49NIdwuvB2xCugt7M9LjGyaT7tpdA/gt6f0OY1n5vRPjgTNz9x3A9Nx9Iemx11WXSx7XzcAGufxE2t8NMa8wzf8FnF8l7iPzPB6S14nnyHdDA2eTHpy30rytNdw8rB8U6k/N030Q8Dzw1jwPZgI7Foa1SSGetYBbgA/UMW0nFL63EFinuD50mM73AL+qY3ntBbxEejLCINIjxovr1GwK65X/Vv5bE+tPPi/pQ7l7DGkjOD3vCe5B+jG8CbgVOJq08b5T6QBiPdofkvYK6YFqFe+RdAJpI78x8KCkW4AhEfHXXOf/SA/Og7QR3aGwh78R6Vk+jxSG+U7gysjPS5JU7TlDz5N+3OdL+h1pY1zNW/Le+lBgMCs/f+eqSA/Oe0jS63LZ+4ALK+OOiKfzHu7bgV+q/YVt69QY39UR8SLwYt5j3o2UjPcG7s51BudpfqzDd28BjlVqo3kIGKb0nKY9gc+TnjlUbbnsQUp+t+bytYHbCsOtPPhuJumBfNXcGOl9DkslPUdKUJAS+A41vlPPcH8TESHpftKjx+8HkPQgKXHeU+U73ycl5t9IOqCLabus0H0f8HNJVwFXVRnuSGBJh7Jqy+tZ4I6ImJtj/QVpGV6Rv7OYlGCeqjHNA5oTRz8haS/SBnHPiHhB0k3Aurn3ZaSX0DxM2liH0i/04og4ucrgXoqIV/Jw1yXt9Y6LiPmSTsvD7ex1mCLtqXb1ALVOn3cTESsk7UZ6UdN44Bjg36tUvQj4YETcK+lI0t5kxcsd4qr87zjuNYBnI2LHLmKuFnflkdXfjoifFnsovT+hvWLEAqWG131Je9kbk5bNsohYWmu5SPoA6SU8tU6hVKbzFWr/rovz4tXC51fr+E49wy0Os+Zw8zLagrQ8Ic27zqbtn4Xu95OO5A4EviJp+2h/MRXAi7Sv9xXVlldn5eRhvFgjngHPbRz9x0bAMzlpvJm0h1rxa1KD8mG0773dABwiaVN47f3FW1QZbuVH+GTeKz8EINIbxJbmIxlIG/aKa4HPKj0KGklvVHqybtHNwIckrZfPS3+g44jz+DaKiGuAY0nvagBYSjrlUjEEWJTHd3iVaejoOuATam8L2TjS+w0ekXRoLpOkt9X4/kGS1pU0nJSk7szT/InCuflRed52jBXS3vSxeR7cAhyf/0Pt5TIdeIekrXP5+pLeWMe09iqSdiFN78fykSDUOW25DWJMRNwInED7EWbRLGDrDmXVlhfAbkpP110D+CjwlzweAa8nnSazKpw4+o8/AGtKug/4BunHCLy2kX8I2CIi7shlDwFfJr3Z7z7SOd5VXvEa6UqY80inM66i/UcH6ammkyXdRtprfC6Xn5/Hd5dSg/lP6bDnGelVmJeRTmP8ivYNZ9EQ4Lc5vj8DX8zlU4H/p9SwuxXwFdIb0a4nHVV1KiL+QGoAnSHpHtKGDFLSmSip8iTgWq8ivQP4HWkefyMiFkbEdaTTdbflUzZXkE7lPUU6BfOApDPz928hvfd7DunJqhtXpr/WcomIJaS2hV/k8ulArcb73uwY0vTemBvIzy8xbYOAn+X5ezdwdnS4UisiHgY2yjsjFassr1x+G7ltjnQa9cpcvgup7al4JGMFfjqurTZJgyNiWe4+ibSB+0KLw7IBTtIXgaURcX4ndfYiXYxxQJV+3wemRcQNjYuyb/MRh3XH+/Ne4wOkxu7TWx2QGek1tC93Wau2B5w0OucjDjMzK8VHHGZmVooTh5mZleLEYWZmpThxmJlZKU4cZmZWyv8HwMn2kh0BzlYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "h_L2=plt.hist([x[1] for x in shimmer_count_L2], \n",
    "              bins=200, range=(0, 2000), alpha=0.75)\n",
    "plt.xlabel(\"average distance between minimizers (bp)\")\n",
    "plt.ylabel(\"read count\")\n",
    "plt.legend( (\"level 2\",) )\n",
    "plt.title(\"Distribution of the distance between shimmers\")\n",
    "print(\"mean distance:\", np.mean([x[1] for x in shimmer_count_L2]) )\n",
    "print(\"median distance:\", np.median([x[1] for x in shimmer_count_L2]) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(229.0, 8190)"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k_count = []\n",
    "for rid in all_rids[:50000]:\n",
    "    if rid not in mmer_index_L2:\n",
    "        continue\n",
    "    shimmers_L2, kset_L2 = get_shimmers_for_read(rid, mmer_index_L2, py_mmer_L2)\n",
    "    if len(shimmers_L2) < 2:\n",
    "        continue\n",
    "    k_count.append( (np.mean([c[-1] for c in shimmers_L2]), rid ) )\n",
    "\n",
    "k_count.sort()\n",
    "k_count[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean count: 79.56188576576095\n",
      "median count: 77.73333333333333\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcXFWd9/HP10TAlRAJTEgCCRBF5EGWlkUZXIIIDBIeZYk6EBbNMLI5yAiI84CoI7igoiM8LUQCIhFxISgiEWFAZTFhCxCBDMGkSSSBQEAQMOE3f9xTSaVT1X1vd61d3/frVa+6de659/5uVXf96pxzF0UEZmZmeb2q2QGYmVl7ceIwM7NCnDjMzKwQJw4zMyvEicPMzApx4jAzs0KcOKxfki6S9B81WteWkv4qaVh6fbOkj9di3Wl9v5I0tVbrK7DdL0p6UtJfCi73mKR9qsz7R0kP1SZCs9oZ3uwArLkkPQZsDqwCVgMPApcB3RHxCkBEHFdgXR+PiN9UqxMRi4DXDy7qNds7G9g2Iv65bP3712LdBeMYB3wa2CoiltVqvRFxK/CWWq1vqJE0HlgIvDoiVjU3ms7iFocBfDAi3gBsBZwLnAZcUuuNSBqqP1S2Ap6qZdJoZUP4c7ScnDhsjYhYGRGzgMOBqZJ2AJB0qaQvpulNJf1C0jOSVki6VdKrJF0ObAlcm7qiPiNpvKSQdKykRcBvy8rKv3y2kXSnpJWSrpE0Mm3rPZJ6ymMsde1I2g/4LHB42t69af6arq8U1+ck/VnSMkmXSdo4zSvFMVXSotTNdGa190bSxmn55Wl9n0vr3weYDWyR4ri0wrIV37OyKjtJui/t/48kbVRp/9O+/3uq+7ykSyRtnrrnnpP0G0mb9Nq/oyUtlvS0pOMkvSMt/4yk7/SK8xhJ81PdX0vaqmxeSDpe0iPAI1Xeo70k/SGte7Gko/p679K8syX9oGwd6/x9pM/zC5J+n/bxBkmbpuq3pOdn0nu/p6RtJf13ei+flPSjap+pDZwTh60nIu4EeoB/rDD702neKLIurs9mi8QRwCKy1svrI+IrZcu8G3gr8IEqmzwSOAbYgqzL7IIcMV4P/Cfwo7S9t1eodlR6vBfYmqyL7Du96uxF1h00Cfh/kt5aZZPfBjZO63l3ivno1C23P7AkxXFUhWUrvmdl8w8D9gMmADummKv5MPB+4M3AB4FfpfVtSvb/fFKv+rsDE8l+DHwTOBPYB3gbcJikdwNIOjit50MpzluBK3ut6+C0vu17ByVpyxTLt9PyOwH3pNkV37s+9rG3j6b6mwEbAKem8r3T84j03t8GfAG4AdgEGJu2bTXmxGHVLAFGVij/OzCarD//7xFxa/R/wbOzI+L5iPhblfmXR8T9EfE88B9kX2jDBh76Gh8Dzo+IRyPir8AZwJRerZ3PR8TfIuJe4F5gvQSUYjkcOCMinouIx4CvA0fkjKO/9+yCiFgSESuAa8m+dKv5dkQ8ERGPk3253xERd0fES8DPgJ171f9CRLwYETcAzwNXRsSysuVL9f8F+HJEzE/jBf9J1hLaqmxdX46IFVU+x48Bv4mIK9M+PhUR99TgvQP4fkQ8nLZ7FX2/P38n6zrcIu337wpsx3Jy4rBqxgArKpR/FVgA3CDpUUmn51jX4gLz/wy8muwX9GBtkdZXvu7hZL/6S8qPgnqBygP3m5L90u29rjE54+jvPcsTQ8kTZdN/q/C697J5628FfCt1Mz1D9tmLdfexr89xHPA/FcoH+95BsffnM2Rx3ynpAUnHFNiO5eTEYeuR9A6yf+z1fq2lX42fjoitybpKTpE0qTS7yir7a5GMK5vekuxX45Nkv5BfWxbXMLJukLzrXUL2hVi+7lWs++WZx5Os/SVbvq7H8yzcz3vWKhYD/xIRI8oer4mIP5TV6ev9XgxsU6G8v/dunc8Y+IcCMa8XT0T8JSI+ERFbkLWivitp2wLrtBycOGwNSW+UdCAwE/hBRMyrUOfANAAp4FmyQ3hXp9lPkPVjF/XPkraX9FrgHODqiFgNPAxsJOmfJL0a+BywYdlyTwDjew00l7sS+DdJEyS9nrVjIoUO3UyxXAV8SdIbUvfNKcAP+l4y08971iouAs6Q9DZYM6B9aIHlrwD2kXSYpOGS3iRppxzv3T3A3srO79mYrDsxr+XAK5T9zUk6VNLY9PJpsuTSau9123PiMMiOhHqO7FfjmcD5VB+8nAj8BvgrcBvw3Yi4Oc37MvC51N1xapXlK7kcuJSsS2Ij0gBvRKwEPglcTPYL9XmyQeaSH6fnpyTdVWG909O6byE73v9F4MQCcZU7MW3/UbKW2A/T+vPo6z1rCRHxM+A8YKakZ4H7yQb98y6/CDiA7ECAFWQJoTReVPW9i4jZwI+A+4C5wC8KbPMF4EvA79Pf3B7AO4A7JP0VmAWcHBEL867T8pFv5GRmZkW4xWFmZoXULXFImq7spKv7K8w7NZ3ks2l6LUkXSFqQTk7apazuVEmPpEfDr0FkZmbrqmeL41Kyk5rWoey6Pu8nO1msZH+yfuCJwDTgwlR3JHAW2UlHuwFnKZ0Za2ZmzVG3xBERt1D5PIBvkB1rXT64Mhm4LDK3AyMkjSY703h2OunoabJLO6yXjMzMrHEaerEySQcBj0fEvdmRiWuMYd2Ti3pSWbXySuueRtZa4XWve92u2223XQ0jNzMb+ubOnftkRIzqr17DEkc6Rv9MYN9KsyuURR/l6xdGdAPdAF1dXTFnzpwBRmpm1pkk/bn/Wo09qmobsou43avsvg1jgbsk/QNZS6L87OGxZGf9Vis3M7MmaVjiiIh5EbFZRIyPiPFkSWGXiPgL2Yk6R6ajq/YAVkbEUuDXwL6SNkmD4vumMjMza5J6Ho57JdlZsm+R1CPp2D6qX0d2VukC4HtkZwuTrhb6BeCP6XFOKjMzsyYZkmeOe4zDzKw4SXMjoqu/ej5z3MzMCnHiMDOzQpw4zMysECcOMzMrxInDzMwKceKwppnSfVuzQzCzAXDisKZzAjFrL04c1hS9k8WU7tucQMzahBOHmZkV4sRhZmaFOHGYmVkhDb2Rk5nHMczan1scZmZWiBOHmZkV4sRhZmaFOHFYy/E4iFlrc+IwM7NCnDjMzKwQJw4zMyvEicOayuMZZu3HicPMzApx4rCGyNuycAvErPXVLXFImi5pmaT7y8q+KulPku6T9DNJI8rmnSFpgaSHJH2grHy/VLZA0un1itfMzPKpZ4vjUmC/XmWzgR0iYkfgYeAMAEnbA1OAt6VlvitpmKRhwH8B+wPbAx9Jdc3MrEnqljgi4hZgRa+yGyJiVXp5OzA2TU8GZkbESxGxEFgA7JYeCyLi0Yh4GZiZ6loHcLeVWWtq5hjHMcCv0vQYYHHZvJ5UVq18PZKmSZojac7y5cvrEK4NVCkBDCQR+M6AZq2nKYlD0pnAKuCKUlGFatFH+fqFEd0R0RURXaNGjapNoNY0ThZmravh9+OQNBU4EJgUEaUk0AOMK6s2FliSpquVm5lZEzS0xSFpP+A04KCIeKFs1ixgiqQNJU0AJgJ3An8EJkqaIGkDsgH0WY2M2czM1lW3FoekK4H3AJtK6gHOIjuKakNgtiSA2yPiuIh4QNJVwINkXVjHR8TqtJ4TgF8Dw4DpEfFAvWI2M7P+1S1xRMRHKhRf0kf9LwFfqlB+HXBdDUMzM7NB8JnjZmZWiBOHmZkV4sRhZmaFOHGYmVkhThxmZlaIE4eZmRXixGFtwZcgMWsdThxmZlaIE4e1Dbc6zFqDE4eZmRXixGFmZoU4cZiZWSENvx+HdQ6PSZgNTW5xmJlZIU4cZmZWiBOHmZkV4sRhZmaFeHDcas6D4mZDm1sc1laclMyaz4nDzMwKceIwM7NC+kwckl4l6bBGBWNmZq2vz8QREa8AJwxkxZKmS1om6f6yspGSZkt6JD1vksol6QJJCyTdJ2mXsmWmpvqPSJo6kFjMzKx28nRVzZZ0qqRx6Yt/pKSROZa7FNivV9npwI0RMRG4Mb0G2B+YmB7TgAshSzTAWcDuwG7AWaVkY2ZmzZHncNxj0vPxZWUBbN3XQhFxi6TxvYonA+9J0zOAm4HTUvllERHA7ZJGSBqd6s6OiBUAkmaTJaMrc8RtZmZ10G/iiIgJNdze5hGxNK13qaTNUvkYYHFZvZ5UVq18PZKmkbVW2HLLLWsYsrWa0iG5M6ft2eRIzDpTv11Vkl4r6XOSutPriZIOrHEcqlAWfZSvXxjRHRFdEdE1atSomgZnZmZr5Rnj+D7wMvDO9LoH+OIAt/dE6oIiPS8rW+e4snpjgSV9lJuZWZPkSRzbRMRXgL8DRMTfqNwSyGMWUDoyaipwTVn5kenoqj2AlalL69fAvpI2SYPi+6YyMzNrkjyD4y9Leg2pi0jSNsBL/S0k6Uqywe1NJfWQHR11LnCVpGOBRcChqfp1wAHAAuAF4GiAiFgh6QvAH1O9c0oD5WZm1hx5EsfZwPXAOElXAO8CjupvoYj4SJVZkyrUDdY9aqt83nRgeo44zcysAfIcVXWDpLnAHmRdVCdHxJN1j8zMzFpSnqOqbgR2j4hfRsQvIuLJ0hFWZs3kK+WaNUeewfEJwGmSzior66pTPNbm/GVuNvTlSRzPkI1LbC7pWkkb1zkmMzNrYXkShyJiVUR8EvgJ8Dtgs36WMTOzISrPUVUXlSYi4lJJ86hyBJSZmQ19VROHpDdGxLPAj3tdDXchcGrdI7O24/ENs87QV4vjh8CBwFzWv25Uv1fHNTOzoalq4oiIA9NzLa+Oa2ZmbS7PGAeSdgTGl9ePiJ/WKSYzM2th/SYOSdOBHYEHgFdScQBOHGZmHShPi2OPiNi+7pGYDcCU7tt8QyezBstzHsdtkpw4zMwMyNfimEGWPP5Cdjl1kV3Qdse6RmZmZi0pT+KYDhwBzGPtGIeZmXWoPIljUUTMqnskZmbWFvKMcfxJ0g8lfUTSh0qPukdmlpPPWDdrrDwtjteQjW3sW1bmw3HNzDpUnjsAHt2IQMzMrD30dZHDz0TEVyR9m6yFsY6IOKmukZmZWUvqq8UxPz3PaUQgZmbWHvq6yOG16XlG48IxM7NWl+daVW8mu//GeNa9yOH7BrpRSf8GfJysC2wecDQwGpgJjATuAo6IiJclbQhcBuwKPAUcHhGPDXTbVns+qsmss+Q5qurHZHcBvBhYPdgNShoDnARsHxF/k3QVMAU4APhGRMyUdBFwLHBhen46IraVNAU4Dzh8sHGYmdnA5DmPY1VEXBgRd0bE3NJjkNsdDrxG0nDgtcBS4H3A1Wn+DODgND05vSbNnySp/KZSZmbWQHkSx7WSPilptKSRpcdANxgRjwNfAxaRJYyVZHcZfCYiVqVqPcCYND0GWJyWXZXqv6n3eiVNkzRH0pzly5cPNDwzM+tHnq6qqen538vKBnzrWEmbkLUiJgDPkHWF7V+haukQ4Eqti0qHB3cD3QBdXV3rzTczs9rIcwJgrW8duw+wMCKWA0j6KfBOYISk4alVMRZYkur3AOOAntS1tTGwosYxmZlZTlW7qiS9Lz1/qNJjENtcBOwh6bVprGIS8CBwE3BIqjMVuCZNz2Jtq+cQ4LcR4RaFrcNHdpk1Tl8tjncDvwU+WGHegK9VFRF3SLqa7JDbVcDdZF1MvwRmSvpiKrskLXIJcLmkBWQtjSkD2a6ZmdVGXycAnpWea36tqrTus3oVPwrsVqHui8ChtY7BasO/9M06T54TAEcAR7L+CYC+VpWZWQfKc1TVdcDt+A6A1uJKrZ+Z0/ZsciRmQ1uexLFRRJxS90jMzKwt5DkB8HJJn6jVCYBmZtbe8rQ4Xga+CpzJ2hPvBnwCoA0dHhg360x5EscpwLYR8WS9gzEzs9aXp6vqAeCFegdiZmbtIU+LYzVwj6SbgJdKhT4c18ysM+VJHD9PDzMzs1wXOfStY83MbI08YxxmZmZrOHGYmVkhfV1W/fL0fHLjwjEzs1bXV4tjV0lbAcdI2qT8rHGfOW5m1rn6Ghy/CLie7Azxuax7C1efOW5m1qGqtjgi4oKIeCswPSK2jogJZQ8nDTOzDtXv4HhE/Kukt0s6IT12bERgZoPh62iZ1U+/iUPSScAVwGbpcYWkE+sdmJmZtaY8Z45/HNg9Ip4HkHQecBvw7XoGZq3Nv+jNOlee8zhEdr2qktWsO1BuZmYdJE+L4/vAHZJ+ll4fDFxSv5DMzKyV5blW1fmSbgb2ImtpHB0Rd9c7MLOBcjeaWX3laXEQEXcBd9Vqo5JGABcDO5CdE3IM8BDwI2A88BhwWEQ8LUnAt4ADyO4LclSKx8zMmqBZ16r6FnB9RGwHvB2YD5wO3BgRE4Eb02uA/YGJ6TENuLDx4Vo7csvDrD4anjgkvRHYmzROEhEvR8QzwGSgdAn3GWRjKaTyyyJzOzBC0ugGh21mZkmfiUPSMEm/qfE2twaWA9+XdLekiyW9Dtg8IpYCpOfNUv0xwOKy5XtSmZmZNUGfiSMiVgMvSNq4htscDuwCXBgROwPPs7ZbqpJKh/7GepWkaZLmSJqzfPny2kRqZmbryTM4/iIwT9Jssi95YFD3HO8BeiLijvT6arLE8YSk0RGxNHVFLSurP65s+bHAkt4rjYhuoBugq6trvcRiZma1kSdx/DI9aiIi/iJpsaS3RMRDwCTgwfSYCpybnq9Ji8wCTpA0E9gdWFnq0rLG84CzmeW657ik1wBbpi/6WjiR7JpXGwCPAkeTdZtdJelYYBFwaKp7HdmhuAvIDsc9ukYxmJnZAPSbOCR9EPgasAEwQdJOwDkRcdBANxoR9wBdFWZNqlA3gOMHui0zM6utPIfjng3sBjwDa770J9QxJjMza2F5EseqiFjZq8yDz9ZWPDZjVjt5Bsfvl/RRYJikicBJwB/qG5a1onb88m3HmM1aXZ4Wx4nA24CXgCuBZ4FP1TMoMzNrXXmOqnoBODPdwCki4rn6h2VmZq0qz61j3yFpHnAf2YmA90ratf6hmZlZK8ozxnEJ8MmIuBVA0l5kN3fasZ6BmZlZa8ozxvFcKWkARMTvAHdXmZl1qKotDkm7pMk7Jf1/soHxAA4Hbq5/aGZm1or66qr6eq/XZ5VN+zyODuPDWs2spGriiIj3NjIQMzNrD3muVTUCOJLsXuBr6g/isupmZtbG8hxVdR1wOzAPeKW+4ZiZWavLkzg2iohT6h6JmZm1hTyH414u6ROSRksaWXrUPTKzGvMAv1lt5GlxvAx8FTiTtUdTBbB1vYIyM7PWlSdxnAJsGxFP1jsYMzNrfXm6qh4gu2WrmZlZrhbHauAeSTeRXVod8OG4ncRjA2ZWLk/i+Hl6mJmZ5bofx4xGBGJmZu0hz/04Fkp6tPejEcGZ1Zq73cwGL09XVVfZ9EbAocCgz+OQNAyYAzweEQdKmgDMTOu+CzgiIl6WtCFwGbAr8BRweEQ8Ntjtm5nZwPTb4oiIp8oej0fEN4H31WDbJwPzy16fB3wjIiYCTwPHpvJjgacjYlvgG6memZk1SZ6uql3KHl2SjgPeMJiNShoL/BNwcXotsmR0daoyAzg4TU9Or0nzJ6X6ZmbWBHm6qsrvy7EKeAw4bJDb/SbwGdYmoDcBz0TEqvS6BxiTpscAiwEiYpWklam+T0i0AZnSfRszp+3Z7DDM2laeo6pqel8OSQcCyyJirqT3lIorbTrHvPL1TgOmAWy55ZY1iNTAg8lmtr489+PYEPgw69+P45wBbvNdwEGSDiAbbH8jWQtkhKThqdUxFliS6vcA44AeScOBjYEVvVcaEd1AN0BXV5fvUGhmVid5LjlyDdk4wyrg+bLHgETEGRExNiLGA1OA30bEx4CbgENStalpuwCz0mvS/N9GhBODDcqU7tvcmjIboDxjHGMjYr+6RwKnATMlfRG4G7gklV9Cdmn3BWQtjSkNiMXMzKrIkzj+IOn/RMS8Wm88Im4Gbk7TjwK7VajzItm5I2Zm1gLyJI69gKMkLSS7yKGAiIgd6xqZmZm1pDyJY/+6R2FmZm0jz+G4f25EIGZm1h7yHFVlNmT5yCqz4pw4zMyskDxjHNZh/CvczPriFoeZmRXixGFmZoU4cZiZWSFOHGZmVogTh5mZFeLEYWZmhThxWMfz4cdmxThxmJlZIU4ctg7/+jaz/jhxmJlZIU4cZrilZVaEE4dZ4vuQm+XjxGFmZoU4cdga/rWd8ftg1jcnDjMzK8SJw8zMCnHiMDOzQhqeOCSNk3STpPmSHpB0ciofKWm2pEfS8yapXJIukLRA0n2Sdml0zGZmtlYzWhyrgE9HxFuBPYDjJW0PnA7cGBETgRvTa4D9gYnpMQ24sPEhm5lZScMTR0QsjYi70vRzwHxgDDAZmJGqzQAOTtOTgcsiczswQtLoBodtHchHV5lV1tQxDknjgZ2BO4DNI2IpZMkF2CxVGwMsLlusJ5X1Xtc0SXMkzVm+fHk9wzYz62hNSxySXg/8BPhURDzbV9UKZbFeQUR3RHRFRNeoUaNqFaaZmfXSlMQh6dVkSeOKiPhpKn6i1AWVnpel8h5gXNniY4EljYp1qPNlNirze2JWXTOOqhJwCTA/Is4vmzULmJqmpwLXlJUfmY6u2gNYWerSMjOzxhvehG2+CzgCmCfpnlT2WeBc4CpJxwKLgEPTvOuAA4AFwAvA0Y0N18zMyjU8cUTE76g8bgEwqUL9AI6va1Dmrhkzy81njpv1oTyhOrmaZZw4zMyskGaMcViL8C/ofPw+ma3LLQ6zApxEzJw4zMysICcOswFwy8M6mRNHh/IXn5kNlBOHWUFOutbpnDjMzKwQJw6zGnArxDqJz+MwGyAnC+tUbnF0IH/hmdlgOHGY1YgTsnUKJ44O4ps21V/5e+z32oYqJw4zMyvEiaND+NdvY/ly7DaUKbtP0tDS1dUVc+bMaXYYLcFfWq1l5rQ9mx2CWVWS5kZEV3/13OIwM7NCnDiGMLc2WpMPUrB258Rh1kCVEobHQ6zd+MzxIchfPu0hT8Iolc+ctuc602bN5MHxIWRK923rfMHY0OYEYrWWd3C8bVockvYDvgUMAy6OiHObHFJL6P0r1Emjc1T7rCu1TnrXLdXpqyXjFo5V0xYtDknDgIeB9wM9wB+Bj0TEg5Xqd0KLwwnCBqNay7SvHyDVEkgpAfWetvaTt8XRLoljT+DsiPhAen0GQER8uVL9eiWO8l9olX6dlZeV/+P1909VqR/brNP19z/WX1mlFlPvJFdtW73XUyQxtnNLbagljkOA/SLi4+n1EcDuEXFCWZ1pwLT08i3AQ4PY5KbAk4NYvh112j532v6C97lTDGaft4qIUf1VapcxDlUoWyfjRUQ30F2TjUlz8mTdoaTT9rnT9he8z52iEfvcLudx9ADjyl6PBZY0KRYzs47WLonjj8BESRMkbQBMAWY1OSYzs47UFl1VEbFK0gnAr8kOx50eEQ/UcZM16fJqM522z522v+B97hR13+e2GBw3M7PW0S5dVWZm1iKcOMzMrBAnjjKS9pP0kKQFkk5vdjz1IukxSfMk3SNpTiobKWm2pEfS8ybNjnMwJE2XtEzS/WVlFfdRmQvS536fpF2aF/nAVdnnsyU9nj7reyQdUDbvjLTPD0n6QHOiHjhJ4yTdJGm+pAcknZzKh+zn3Mc+N/Zzjgg/snGeYcD/AFsDGwD3Ats3O6467etjwKa9yr4CnJ6mTwfOa3acg9zHvYFdgPv720fgAOBXZOcL7QHc0ez4a7jPZwOnVqi7ffob3xCYkP72hzV7Hwru72hglzT9BrLLEm0/lD/nPva5oZ+zWxxr7QYsiIhHI+JlYCYwuckxNdJkYEaangEc3MRYBi0ibgFW9Cquto+TgcsiczswQtLoxkRaO1X2uZrJwMyIeCkiFgILyP4H2kZELI2Iu9L0c8B8YAxD+HPuY5+rqcvn7MSx1hhgcdnrHvr+QNpZADdImpsu1QKweUQsheyPE9isadHVT7V9HOqf/Qmpa2Z6WRfkkNpnSeOBnYE76JDPudc+QwM/ZyeOtfq9rMkQ8q6I2AXYHzhe0t7NDqjJhvJnfyGwDbATsBT4eiofMvss6fXAT4BPRcSzfVWtUDZU9rmhn7MTx1odc1mTiFiSnpcBPyNruj5Raran52XNi7Buqu3jkP3sI+KJiFgdEa8A32NtN8WQ2GdJryb7Ar0iIn6aiof051xpnxv9OTtxrNURlzWR9DpJbyhNA/sC95Pt69RUbSpwTXMirKtq+zgLODIddbMHsLLU1dHuevXh/1+yzxqyfZ4iaUNJE4CJwJ2Njm8wJAm4BJgfEeeXzRqyn3O1fW7459zsowRa6UF21MXDZEcenNnseOq0j1uTHWVxL/BAaT+BNwE3Ao+k55HNjnWQ+3klWZP972S/uo6tto9kzfn/Sp/7PKCr2fHXcJ8vT/t0X/oSGV1W/8y0zw8B+zc7/gHs715k3S73AfekxwFD+XPuY58b+jn7kiNmZlaIu6rMzKwQJw4zMyvEicPMzApx4jAzs0KcOMzMrBAnDutIkg7q7wrIkraQdHXB9Z4jaZ/BRddYkkZI+mSz47D24cNxzVqMpOERsaqB2xsP/CIidmjUNq29ucVhQ4qk8ZL+JOliSfdLukLSPpJ+n+7PsFuqd5Sk76TpS9N9Gv4g6VFJh5St6/6y+j+XdK2khZJOkHSKpLsl3S5pZNm6DpHUVXZvhHmSIs3fRtL16QKTt0rarmy58yXdBJzXa5+GSfpaWs99kk5M5ZPS9uelC9ttmMofk7Rpmu6SdHOaPjvVuznt50lpE+cC26RYv1q/T8eGiuHNDsCsDrYFDgWmkV1K5qNkZ9weBHyWypeMH53qbEd25m2lLqodyK5GuhHZ5alPi4idJX0DOBL4ZqliRMwhu+Ac6cv4+jSrGzguIh6RtDvwXeB9ad6bgX0iYnWv7U4ju5fCzhGxStmNijYCLgUmRcTDki4D/rU8hiq2A95Ldi+HhyRdSHbPih0iYqd+ljUDnDhsaFoYEfMAJD0A3BgRIWkeML7KMj+P7AJxD0ravEqdmyK7B8JzklYC16byecCOlRaQdBjZzZX2TVc0fSfw4+ySQ0B2g52SH1dIGgD7ABeVuq9XG0LoAAABIElEQVQiYoWkt6f9fDjVmQEcT/+J45cR8RLwkqRlQLV9NavKicOGopfKpl8pe/0K1f/my5epdCnqwuuV9Dbg88DeEbFa0quAZ/r4Zf98lXKx/qWwq8UIsIq13dAb9ZpXvg+rK8Vt1h+PcZjVgaSNye4ieWRELAeI7L4JCyUdmuootRz6cwNwnKThabmRwJ+A8ZK2TXWOAP47TT8G7JqmP5xj/c+RdV2Z5eLEYVYfBwNbAd8rDZKn8o8Bx0oqXZ04z+2JLwYWAfel5T4aES8CR5N1e80ja/VclOp/HviWpFvJWhV9ioingN+ngwk8OG798uG4ZmZWiFscZmZWiBOHmZkV4sRhZmaFOHGYmVkhThxmZlaIE4eZmRXixGFmZoX8Ly9uOfNY70TFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist( [c[0] for c in k_count], \n",
    "         bins=250, range=(0,250), alpha=0.75 );\n",
    "plt.xlabel(\"minimizer count\");\n",
    "plt.ylabel(\"number of minimizer\");\n",
    "plt.title(\"Distribution of shimmer counts\")\n",
    "print(\"mean count:\", np.mean([c[0] for c in k_count]) )\n",
    "print(\"median count:\", np.median([c[0] for c in k_count]) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "nread_ovlp = []\n",
    "for rid in all_rids[:50000]:\n",
    "    if rid not in mmer_index_L2:\n",
    "        continue\n",
    "    shimmers_L2, kset_L2 = get_shimmers_for_read(rid, mmer_index_L2, py_mmer_L2, mc_h=240)\n",
    "    if len(shimmers_L2) < 2:\n",
    "        continue\n",
    "    rids = set()\n",
    "    n0 = shimmers_L2[0]\n",
    "    for n1 in shimmers_L2[1:]:\n",
    "        v = shimmer_ffi.new(\"mp256_v *\")\n",
    "        shimmer4py.get_shimmer_hits(v, py_mmer_L2, n0[4], 16)\n",
    "        for i in range(v.n):\n",
    "            if  (v.a[i].x1 >> 8) == n1[4]:\n",
    "                rids.add(v.a[i].y0 >> 32)\n",
    "        n0 = n1\n",
    "    nread_ovlp.append(len(rids))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean count: 52.091704599779675\n",
      "median count: 50.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcHFW5//HP14CEPYkJXMhCWIIKiCwji6A3XhUBUUBFwlVJcIkoKCr+BMQriBfFC4LiAgZZFQJBRKKiEBDEJQEmEJKwSdhkSEzCYgARNOH5/XFOJ5VO90z3zPR0T+b7fr36NVWntqerq+qZOqcWRQRmZmb1eFWzAzAzs/7HycPMzOrm5GFmZnVz8jAzs7o5eZiZWd2cPMzMrG5OHp2QdL6k/+mleY2R9IKkQbn/Vkkf74155/n9RtLE3ppfHcv9X0lPSfpbXy+7LI7xkjqauPxDJT2Rf+NdmxhHU9dDb5AUkrbL3Z3ug8VxrTJJY/N6Wqc359urM+tPJD0GbA4sB1YA9wGXAVMi4hWAiDi6jnl9PCJuqjZORPwV2KhnUa9c3qnAdhHx4cL8D+iNedcZx2jgeGCriFjS18tvMWcBx0bEdc0OZG1S6z7YFUljgUeBdSNieW/Mc6Ab6Gce74mIjYGtgDOAE4ALe3shvZ3xW8hWwNNrW+Lo5u+1FXBvb8dSj7V4O1tr1PIb9ZvfMSIG5Ad4DHhHWdkewCvATrn/EuB/c/dw4FfA34FngD+Qku9P8jT/BF4AvgSMBQL4GPBX4LZC2Tp5frcC3wTuAJYB1wHD8rDxQEeleIH9gX8B/87Lu6cwv4/n7lcBXwEeB5aQzqg2zcNKcUzMsT0FnNzJeto0T780z+8ref7vyN/5lRzHJRWmHQ90kM5OlgCLgKMKw1fGnPsnAX8s9AfwaeAh4Hng68C2wEzgOWAa8OqyZX05f6fHgA8V5rUe6ezgr8Bi4Hxg/bJpTwD+BvykwnepuE7zfF/Isf4DeLjKenwzcGf+re8E3pzLJwDtZeN+Hpjenbgp23aAE4GH8/q7Dzi0bH3/CfhejusB4O2dbAujgZ/nbeFp4Pu5fFvgd7nsKeByYEjZtvtFYG5ezlXA4MLw/5e3jYXAR/O63K58H6xh3HcDd+dt4wng1MJ0f83jvpA/ewPbAb/PMT0FXFXle4/N007Oy10EHF+2bZTW89Ok7XJY2bQrjwWd7CerbX/AQcAc0jHnz8DONf6ug/I28xTwCHAMhWNPrx1DG3mAbuUPFZJHYSP7VPmGSzrQnw+smz9vAVRpXoUN5jJgQ2B9KiePJ4Gd8jjXAD8tbkzV4gVOLY1bGH4rq5LHR4EFwDakqrKfFzbIUhwX5LjeCLwMvL7KerqMlNg2ztP+BfhYtTgr7BTLgdPyOjsQeBEYWh5z7p/EmsljOrAJsGOO8+b8vTbNO83EsmWdTTrg/ifpYP7aPPw7eV7D8nf5JfDNsmm/laddv8J3qbpOC7FuV2U9DAOeBT5Cqio+Ive/BtiAdAAYVxj/TmBCd+Iu/02Aw4AtSQe4w/M62aKwvpeTktW6efgy8oGv7DsMAu4BziFtr4OBffOw7YB35hhGkP5Z+k7ZtntHjmMYcD9wdB62PykplvaDK6iSPGoYdzzwhvxdd87jHlK23a9TiGsqcHIef+X3qfDdS9NOzct9AymBlvbHzwGzgFF5HfwImFrtWNDJflL8HXcj/ZOyZ173E/N6XK+G3/Vo0j8Co/P6vqX8u/fKMbSRB+hW/lA9ecwi/ydetuGeRjqIrnGAKJ9XYYPZpkJZMXmcURi+A+mMYhA9Tx43A58uDHst6UxlnUIcowrD7yAfrMrmOYh0wN6hUPZJ4NbCRt9V8vgnq++wS4C9ymPO/ZNYM3nsU+ifDZxQ6P82+SDFqh1ww8LwacD/AMo717aFYXsDjxam/ReF/4YrfJeq67QQa7Xk8RHgjrKymcCk3P1T4Ku5exwpmWzQnbhr+E3mAAcX1vdC8j9BhW3hIxWm25t0wOzyAAQcAtxdtu1+uND/f8D5ufsiVt8Ptqd68uh03ApxfAc4p9L+l8suA6ZQ2BeqzKc07evKvsOFuft+CmdswBasub9t08n8K/2O5wFfLxvvQeA/a/hdf0dOzrl/v/Lv3hufgd7mUclIUrVUuTNJ/3neKOkRSSfWMK8n6hj+OOm/v+E1Rdm5LfP8ivNeh3SBQEnx6qgXqdyYPxx4dYV5jawjlqdj9QbKasuqZnGh+58V+ovzejYi/lHof5y0LkaQDsazJf1d0t+B3+bykqUR8VIncdSyTmudtjR9aT1eQTobAfhv4BcR8WJvxC3pSElzCtPvxOrb2JORjzCFuLasMKvRwONRobFZ0maSrpT0pKTnSMmwfDuutr1tyZr7QTWdjitpT0m3SFoqaRnpP/DO9qcvkRL0HZLulfTRTsalwrJL62kr4NrCOr6fdBHO5lWmraT8d9wKOL40zzzf0aVldvG71rNOu83Jo0DSm0g79B/Lh0XE8xFxfERsA7wH+IKkt5cGV5lltfKS0YXuMaT/Vp4i/be5QSGuQax+wOhqvgtJG19x3stZ/cBbi6dyTOXzerLO+VSz2vcE/qOH8xsqacNC/xjSuniKlGh2jIgh+bNpRBQTTyPXafm0pelL6/FGYLikXUhJ5Ipc3qO4JW1Fqp48FnhNRAwB5pMOmCUjJRX7S+us3BPAmCqNud/McewcEZsAHy5bRmcWseZ+0N1xryBV8Y2OiE1J1cylONZYTxHxt4j4RERsSTqj/mEXl/2WL7u0np4ADij8RkMiYnBEFPeTrrav8uFPAKeXzXODiJhaw+9azzrtNicPQNImkg4CriRVB82rMM5BkrbLO9pzpP8sVuTBi0l14fX6sKQdJG1Aqhb7WUSsILUrDJb0bknrkhpq1ytMtxgYK6na7zcV+LykrSVtBHyD1BhY1yWKOZZpwOmSNs4b7RdI/1n2hjnA+yRtkHfaj/XCPL8m6dWS3kJqcLw60qXXFwDnSNoMQNJISe+qY749WafXA9tL+m9J60g6nFRN+SuAPI+fkc5uhwEzcnlP496QdFBamqc9ivQfatFmwGclrSvpMOD1Od5yd5AOSmdI2lDSYEn75GEbkxqh/y5pJKlRu1bTgEmF/eCUHoy7MfBMRLwkaQ/SWVzJUtLFHSv3U0mHSRqVe58lrasVVPc/eVvdETiK1PAPKUmdnvcPJI2QdHAn86nFBcDR+WxKeZ2/W9LGdP27TiP9pqMkDSU1rve6gZ48finpeVKWP5nU2HpUlXHHATeRdpKZwA8j4tY87JvAV/Ip5BfrWP5PSHW6fyM12H0WICKWka4y+jHpv9N/kK7GKLk6/31a0l0V5ntRnvdtpGvbXwI+U0dcRZ/Jy3+EdEZ2RZ5/bziHVNe7GLiUdJVOT/yNdBBYmOd1dEQ8kIedQKp2nJWrVm4itVvUqtvrNCKeJiWy40lX43wJOCginiqMdgXpCraryxJSt+OOiPtI7UIzSev4DaSrq4puJ23bTwGnAx/I8ZbPawXpjHs70kUlHaSGWoCvkRp4lwG/Jl1MUJOI+A2pbeJ3+Xv+rgfjfho4Le/TXyUdREvTvpi/35/yfroX8CbgdkkvkM5YjouIRzsJ9/d5uTcDZ0XEjbn8u3n6G/OyZ5EaurstItqBTwDfJ23TC0htVLX8rhcAN5AucLiLOn6PepSuFjKzAUbSJNIFC/s2O5ZW5hsMKxvoZx5mZtYNTh5mZla3hiUPSaPzZXP358vgjsvlwyTNkPRQ/js0l0vSuZIWSJorabfCvCbm8R9SEx7+Z7Y2iohLXGXVtYh4LCLkKqvVNazNQ9IWpDse78pXCMwm3Tw0iXRFxBlK90oMjYgTJB1IaoA8kNTY9N2I2FPSMKAdaCNdYTAb2D0inm1I4GZm1qWGPYArIhaRLu0jIp6XdD/pHoqDSXdUQrrC5lbSFSUHA5flG5ZmSRqSE9B4YEZEPAMgaQbpMQVTO1v+8OHDY+zYsb37pczM1mKzZ89+KiJGdD1mHz2SPV+tsCvpssDNc2IhIhaVrl8nJZbiXZEduaxaeafGjh1Le3t7j2M3MxsoJNV8N3rDG8zzDVXXAJ+LiOc6G7VCWXRSXmlZkyW1S2pfunRp/cGamVlNGpo88t3R1wCXR0TpRpXFuTqq1C5SehdEB6vfUj+KdLNXtfI1RMSUiGiLiLYRI2o68zIzs25o5NVWIr1Y6f6IOLswaDrp8cLkv9cVyo/MV13tBSzL1Vs3APtJGpqvzNovl5mZWZM0ss1jH9KjqOdJmpPLvkx6Y980SaWXoxyWh11PutJqAempm0cBRMQzkr5OescBwGmlxnMzM2uOtfbxJG1tbeEGczOz2kmaHRFttYzrO8zNzKxuTh5mZlY3Jw8zM6ubk4eZmdXNycP6lQlTZjJhysxmh2E24PXJ40nMesoJw6y1+MzDzMzq5jMPa2k+4zBrTU4e1pKcNMxam6utzMysbk4e1i/5qiuz5nLyMDOzurnNw1qKzybM+gefeZiZWd2cPMzMrG5OHtavueHcrDmcPMzMrG5OHmZmVreGJQ9JF0laIml+oewqSXPy57HSu80ljZX0z8Kw8wvT7C5pnqQFks6VpEbFbGZmtWnkpbqXAN8HLisVRMThpW5J3waWFcZ/OCJ2qTCf84DJwCzgemB/4DcNiNfMzGrUsDOPiLgNeKbSsHz28EFgamfzkLQFsElEzIyIICWiQ3o7VjMzq0+z2jzeAiyOiIcKZVtLulvS7yW9JZeNBDoK43TkMjMza6Jm3WF+BKufdSwCxkTE05J2B34haUegUvtGVJuppMmkKi7GjBnTi+Fao5Qus71y8t5NjsTM6tHnZx6S1gHeB1xVKouIlyPi6dw9G3gY2J50pjGqMPkoYGG1eUfElIhoi4i2ESNGNCJ8MzOjOdVW7wAeiIiV1VGSRkgalLu3AcYBj0TEIuB5SXvldpIjgeuaELOZmRU0rNpK0lRgPDBcUgdwSkRcCExgzYbytwKnSVoOrACOjohSY/unSFdurU+6yspXWq2FfJe4Wf/SsOQREUdUKZ9Uoewa4Joq47cDO/VqcGZm1iO+w9zWCn7GlVnfcvKwtYqTiFnfcPIwM7O6OXmYmVndnDzMzKxuTh5mZlY3Jw8zM6ubk4eZmdXNycPMzOrm5GFmZnVz8jAzs7o5eZiZWd2a9TIoG+D8CBGz/s1nHmZmVjcnDzMzq5uTh5mZ1c3Jw8zM6ubkYWZmdWtY8pB0kaQlkuYXyk6V9KSkOflzYGHYSZIWSHpQ0rsK5fvnsgWSTmxUvGZmVrtGnnlcAuxfofyciNglf64HkLQDMAHYMU/zQ0mDJA0CfgAcAOwAHJHHNTOzJmrYfR4RcZuksTWOfjBwZUS8DDwqaQGwRx62ICIeAZB0ZR73vl4O18zM6tCMNo9jJc3N1VpDc9lI4InCOB25rFq5mZk1UV8nj/OAbYFdgEXAt3O5KowbnZRXJGmypHZJ7UuXLu1prNYAE6bM9N3lZmuBPk0eEbE4IlZExCvABayqmuoARhdGHQUs7KS82vynRERbRLSNGDGid4M3M7OV+jR5SNqi0HsoULoSazowQdJ6krYGxgF3AHcC4yRtLenVpEb16X0Zs/VPPsMxa6yGNZhLmgqMB4ZL6gBOAcZL2oVU9fQY8EmAiLhX0jRSQ/hy4JiIWJHncyxwAzAIuCgi7m1UzGZmVptGXm11RIXiCzsZ/3Tg9Arl1wPX92JoZmbWQ77D3MzM6ubkYWZmdXPyMDOzujl5mJlZ3Zw8zMysbk4eZmZWNycPMzOrm5OHmZnVra7kIWmopJ0bFYyZmfUPXSYPSbdK2kTSMOAe4GJJZzc+NDMza1W1nHlsGhHPAe8DLo6I3YF3NDYsMzNrZbUkj3Xy03A/CPyqwfGY9So/XdesMWpJHqeRnmq7ICLulLQN8FBjwzIzs1bW5VN1I+Jq4OpC/yPA+xsZlK19/N+/2dqlavKQ9D06eeVrRHy2IRGZmVnL66zaqh2YDQwGdiNVVT1Eev/4isaHZmZmrarqmUdEXAogaRLwtoj4d+4/H7ixT6IzM7OWVEuD+ZbAxoX+jXKZmZkNULUkjzOAuyVdIukS4C7gG11NJOkiSUskzS+UnSnpAUlzJV0raUguHyvpn5Lm5M/5hWl2lzRP0gJJ50pS3d/SzMx6VZfJIyIuBvYErs2fvUtVWl24BNi/rGwGsFNE7Az8BTipMOzhiNglf44ulJ8HTAbG5U/5PM3MrI/V+myrl4FFwLPA9pLe2tUEEXEb8ExZ2Y0RsTz3zgJGdTaPfHPiJhExMyICuAw4pMaYzcysQbq8z0PSx4HjSAf6OcBewEzgv3q47I8CVxX6t5Z0N/Ac8JWI+AMwEugojNORy8zMrIlqOfM4DngT8HhEvA3YFVjak4VKOhlYDlyeixYBYyJiV+ALwBWSNgEqtW9UvfdE0mRJ7ZLaly7tUYhmZtaJWpLHSxHxEoCk9SLiAeC13V2gpInAQcCHclUUEfFyRDydu2cDDwPbk840ilVbo4CF1eYdEVMioi0i2kaMGNHdEM3MrAu1JI+OfFXUL4AZkq6jkwN4ZyTtD5wAvDciXiyUj5A0KHdvQ2oYfyQiFgHPS9orX2V1JHBdd5ZtZma9p5ZnWx2aO0+VdAuwKfDbrqaTNBUYDwyX1AGcQrq6aj1SEgKYla+seitwmqTlpLvXj46IUmP7p0hXbq0P/CZ/zOpSfLbWlZP3bmIkZmuHLpMHgKR9gXERcbGkEaRG60c7myYijqhQfGGVca8BrqkyrB3YqZY4rfX4gYhma6da3iR4CqmqqXRPxrrATxsZlJmZtbZa2jwOBd4L/AMgIhay+uNKzMxsgKml2upfERGSAkDShg2OydYCrq4yW7vVcuYxTdKPgCGSPgHcBFzQ2LDMzKyV1XK11VmS3km68/u1wFcjYkbDIzMzs5bVafLI917cEBHvID3U0MzMrPNqq4hYAbwoadM+isfMzPqBWhrMXwLmSZpBvuIK/A5zM7OBrJbk8ev8MTMzA2prMK/lxU9mZjaA1PoyKDMzs5WcPGzAmTBlpm9iNOuhqslD0k/y3+P6LhwzM+sPOjvz2F3SVsBHJQ2VNKz46asAzcys9XTWYH4+6b0d2wCzWf2VsJHLzcxsAKp65hER50bE64GLImKbiNi68HHiMDMbwGq5VPdTkt4IvCUX3RYRcxsblpmZtbJaXgb1WeByYLP8uVzSZxodmJmZta5a7jD/OLBnRPwDQNK3gJnA9xoZmJmZta5a7vMQsKLQv4LVG8+rTyhdJGmJpPmFsmGSZkh6KP8dmssl6VxJCyTNlbRbYZqJefyHJE2s7auZmVmj1JI8LgZul3SqpFOBWcCFNc7/EmD/srITgZsjYhxwc+4HOAAYlz+TgfMgJRvgFGBPYA/glFLCMTOz5ugyeUTE2cBRwDPAs8BREfGdWmYeEbfl6YoOBkrPy7oUOKRQflkks0hvLtwCeBcwIyKeiYhnSe8VKU9IZnXzneZm3VdLmwcRcRdwVy8tc/OIWJTnu0jSZrl8JPBEYbyOXFatfA2SJpPOWhgzZkwvhWtmZuVa6dlWldpRopPyNQsjpkREW0S0jRgxoleDMzOzVZqRPBbn6ijy3yW5vAMYXRhvFLCwk3IzM2uSTpOHpEGSburlZU4HSldMTQSuK5Qfma+62gtYlqu3bgD2y8/XGgrsl8vMzKxJOm3ziIgVkl6UtGlELKt35pKmAuOB4ZI6SFdNnQFMk/Qx4K/AYXn064EDgQXAi6RGeiLiGUlfB+7M450WEeWN8NYi3ABtNjA09B3mEXFElUFvrzBuAMdUmc9FwEU1xGpmZn3A7zA3M7O61fQOc0nrA2Mi4sE+iMnMzFpcLQ9GfA8wh/RuDyTtIml6owMzM7PWVculuqeSHgvyd4CImANs3cCYzMysxdWSPJZXuNKq4k16ZmY2MNTSYD5f0n8DgySNAz4L/LmxYZmZWSur5czjM8COwMvAVOA54HONDMrMzFpbLVdbvQicnF8CFRHxfOPDMjOzVlbL1VZvkjQPmEu6WfAeSbs3PjSzvuFHs5vVr5Y2jwuBT0fEHwAk7Ut6QdTOjQzM+hcffM0GllraPJ4vJQ6AiPgj4KorM7MBrOqZR+Ed4ndI+hGpsTyAw4FbGx+amZm1qs6qrb5d1n9Kodv3eZiZDWBVk0dEvK0vAzEzs/6jywZzSUOAI4GxxfFreSS7mZmtnWq52up6YBYwD3ilseGYmVl/UEvyGBwRX2h4JGZm1m/UcqnuTyR9QtIWkoaVPg2PzMzMWlYtyeNfwJnATGB2/rR3d4GSXitpTuHznKTPSTpV0pOF8gML05wkaYGkByW9q7vLNuuM7zQ3q10t1VZfALaLiKd6Y4H5bYS7AEgaBDwJXAscBZwTEWcVx5e0AzCB9HDGLYGbJG0fESt6Ix7rGR9szQamWs487gVebNDy3w48HBGPdzLOwcCVEfFyRDwKLCC9nMrMzJqkljOPFcAcSbeQHssO9NqluhNId66XHCvpSFK12PER8SwwknS1V0lHLluDpMnAZIAxY8b0QnhmZlZJLWcevwBOJ70Aanbh0yOSXg28F7g6F50HbEuq0lrEqjvcVWHyine4R8SUiGiLiLYRI0b0NEQzM6uilvd5XNqgZR8A3BURi/NyFpcGSLoA+FXu7QBGF6YbBSxsUExmZlaDWu4wf5QK/+lHxDY9XPYRFKqsJG0REYty76HA/Nw9HbhC0tmkBvNxwB09XLaZmfVALW0ebYXuwcBhQI/u85C0AfBO4JOF4v+TtAspUT1WGhYR90qaBtwHLAeO8ZVWZmbNVUu11dNlRd+R9Efgq91daH617WvKyj7Syfink9pdzMysBdRSbbVbofdVpDORjRsWkVmTle5duXLy3k2OxKx11VJtVXyvx3JSldIHGxKNmZn1C7VUW/m9HmZmtppaqq3WA97Pmu/zOK1xYZmZWSurpdrqOmAZ6cbAl7sY18zMBoBakseoiNi/4ZGYmVm/UcvjSf4s6Q0Nj8TMzPqNWs489gUm5TvNXyY9ayoiYueGRmZmZi2rluRxQMOjMDOzfqWWS3U7e9eGDVB+CZTZwFZLm4eZmdlqnDzMzKxuTh5mZlY3Jw8zM6ubk4eZmdXNycPMzOrm5GFmZnWr5SZBs5UG0v0dfimUWXVNO/OQ9JikeZLmSGrPZcMkzZD0UP47NJdL0rmSFkiaW/Z2QzMz62PNrrZ6W0TsEhFtuf9E4OaIGAfcnPshPSJlXP5MBs7r80htwJowZeaAOuMyq0Wzk0e5g4FLc/elwCGF8ssimQUMkbRFMwI0M7PmJo8AbpQ0W9LkXLZ5RCwCyH83y+UjgScK03bkstVImiypXVL70qVLGxi6mdnA1swG830iYqGkzYAZkh7oZFxVKIs1CiKmAFMA2tra1hhuZma9o2lnHhGxMP9dAlwL7AEsLlVH5b9L8ugdwOjC5KOAhX0XrZmZFTUleUjaUNLGpW5gP2A+MB2YmEebSHp/Orn8yHzV1V7AslL1lpmZ9b1mVVttDlwrqRTDFRHxW0l3AtMkfQz4K3BYHv964EBgAfAicFTfhzyw+WojMytqSvKIiEeAN1Yofxp4e4XyAI7pg9DMzKwGrXaprpmZ9QNOHmZmVjcnDzMzq5sfjGidckO5mVXiMw8zM6ubk4eZmdXNycPMzOrmNg+zGpW3//glUTaQ+czDzMzq5uRhZmZ1c/IwM7O6OXmYmVnd3GBuFfnmQDPrjM88zMysbk4eZmZWNycPMzOrm5OHmZnVzcnDrJsmTJnpCwtswOrz5CFptKRbJN0v6V5Jx+XyUyU9KWlO/hxYmOYkSQskPSjpXX0d80DiA6KZ1aIZl+ouB46PiLskbQzMljQjDzsnIs4qjixpB2ACsCOwJXCTpO0jYkWfRm1mZiv1+ZlHRCyKiLty9/PA/cDITiY5GLgyIl6OiEeBBcAejY/UzMyqaWqbh6SxwK7A7bnoWElzJV0kaWguGwk8UZisgyrJRtJkSe2S2pcuXdqgqM3MrGnJQ9JGwDXA5yLiOeA8YFtgF2AR8O3SqBUmj0rzjIgpEdEWEW0jRoxoQNRrL7d1dJ/XnQ1ETUkektYlJY7LI+LnABGxOCJWRMQrwAWsqprqAEYXJh8FLOzLeM3MbHV93mAuScCFwP0RcXahfIuIWJR7DwXm5+7pwBWSziY1mI8D7ujDkNdq/o/ZzLqjGVdb7QN8BJgnaU4u+zJwhKRdSFVSjwGfBIiIeyVNA+4jXal1TF9daVU6sPqNcWZmq+vz5BERf6RyO8b1nUxzOnB6w4Iy6wX+Z8MGEt9hbmZmdXPyMDOzuvllUAOUG8rNrCd85jFA+F4EM+tNPvMYYJxAGs8N5zYQ+MzDzMzq5jOPtZzPNJrHZyC2NnPyWAs5YZhZo7naai3gxnAz62s+81iLOIG0Jldf2drIZx5mZlY3n3n0Yz7TMLNmcfLoR5ws1i6uzrL+zMmjBt7JrTc4+dvaxMmjhflgY2atysnDrMlq/SfBZ77WSpw8WojPNKwz5duHk4k1k5NHEzlZWG+oNanU03bndj7rSr9JHpL2B74LDAJ+HBFn9HUMXe2k1Xa48nInDesNXW1H1Yb3ZHusNm1Pk0yt83FSax2KiGbH0CVJg4C/AO8EOoA7gSMi4r5q07S1tUV7e3uPluuDvFnPVEtQtSaursbr7j9kXZ2ddbaMehNXbyXGng6vhaTZEdFW07j9JHnsDZwaEe/K/ScBRMQ3q03Tk+ThpGFm/VVfJY/+Um01Enii0N8B7Fk+kqTJwOTc+4KkB7u5vOHAU92cti+0enzgGHtDq8cHrR9jq8cHvRzjVZ/s0eRb1Tpif0keqlC2xilTREwBpvR4YVJ7rdm3GVo9PnCMvaHV44PWj7HV44P+EWMl/eXBiB3A6EL/KGBhk2IxMxvw+kvyuBMYJ2lrSa8GJgDTmxyTmdmA1S+qrSJiuaRjgRtIl+peFBH3NnCRPa76arBWjw8cY29o9fig9WNs9figf8S4hn5xtZWZmbWW/lJtZWZmLcTJw8zM6ubkUSBpf0kPSlog6cRmxwMgabSkWyTdL+lEyx3HAAAKkklEQVReScfl8mGSZkh6KP8d2uQ4B0m6W9Kvcv/Wkm7P8V2VL3RoZnxDJP1M0gN5Xe7dSutQ0ufz7ztf0lRJg5u9DiVdJGmJpPmFsorrTMm5ed+ZK2m3JsZ4Zv6d50q6VtKQwrCTcowPSnpXs2IsDPuipJA0PPc3ZT12h5NHlh+B8gPgAGAH4AhJOzQ3KgCWA8dHxOuBvYBjclwnAjdHxDjg5tzfTMcB9xf6vwWck+N7FvhYU6Ja5bvAbyPidcAbSbG2xDqUNBL4LNAWETuRLgqZQPPX4SXA/mVl1dbZAcC4/JkMnNfEGGcAO0XEzqTHGp0EkPebCcCOeZof5v2+GTEiaTTpkUt/LRQ3az3WzcljlT2ABRHxSET8C7gSOLjJMRERiyLirtz9POmgN5IU26V5tEuBQ5oTIUgaBbwb+HHuF/BfwM/yKM2ObxPgrcCFABHxr4j4Oy20DklXPq4vaR1gA2ARTV6HEXEb8ExZcbV1djBwWSSzgCGStmhGjBFxY0Qsz72zSPeFlWK8MiJejohHgQWk/b7PY8zOAb7E6jc8N2U9doeTxyqVHoEyskmxVCRpLLArcDuweUQsgpRggM2aFxnfIe0Er+T+1wB/L+zAzV6X2wBLgYtz1dqPJW1Ii6zDiHgSOIv0H+giYBkwm9ZahyXV1lmr7j8fBX6Tu1smRknvBZ6MiHvKBrVMjF1x8lilpkegNIukjYBrgM9FxHPNjqdE0kHAkoiYXSyuMGoz1+U6wG7AeRGxK/APml/Nt1JuNzgY2BrYEtiQVH1RrmW2xwpa7TdH0smkat/LS0UVRuvzGCVtAJwMfLXS4AplLfm7O3ms0rKPQJG0LilxXB4RP8/Fi0uns/nvkiaFtw/wXkmPkar6/ot0JjIkV8FA89dlB9AREbfn/p+RkkmrrMN3AI9GxNKI+Dfwc+DNtNY6LKm2zlpq/5E0ETgI+FCsupmtVWLclvSPwj15vxkF3CXpP2idGLvk5LFKSz4CJbcfXAjcHxFnFwZNBybm7onAdX0dG0BEnBQRoyJiLGmd/S4iPgTcAnyg2fEBRMTfgCckvTYXvR24jxZZh6Tqqr0kbZB/71J8LbMOC6qts+nAkflqob2AZaXqrb6m9OK4E4D3RsSLhUHTgQmS1pO0NalR+o6+ji8i5kXEZhExNu83HcBueTttmfXYpYjwJ3+AA0lXZzwMnNzseHJM+5JOW+cCc/LnQFK7ws3AQ/nvsBaIdTzwq9y9DWnHXABcDazX5Nh2AdrzevwFMLSV1iHwNeABYD7wE2C9Zq9DYCqpDebfpAPcx6qtM1J1yw/yvjOPdOVYs2JcQGo3KO0v5xfGPznH+CBwQLNiLBv+GDC8meuxOx8/nsTMzOrmaiszM6ubk4eZmdXNycPMzOrm5GFmZnVz8jAzs7o5eVi3SbpVUlsfLOez+Um4l3c9dq8sb2ylJ6C2Cknjterpxe9VlSdAS3qhi/kMkfTpRsTYSlr99+yvnDysKQp3Ttfi08CBkW4+bKg642q6iJgeEWd0c/IhpHXbVJ2t83yznI9TLcg/ylou/9d1v6QLlN4XcaOk9fOwlWcOkobnRyUgaZKkX0j6paRHJR0r6Qv5oYKzJA0rLOLDkv6s9B6KPfL0G+Z3GNyZpzm4MN+rJf0SuLFCrF/I85kv6XO57HzSzXLTJX2+bPzBki6WNC8v5225/HZJOxbGu1XS7t2JK6+/P0i6K3/enMvHS7pN6X0R90k6v9JBTtKb8vq5R9IdkjbuYp63atV7Ry7Pd5yX3jXzgKQ/Au8rzH+SpO/n7q0lzczf7+uFcTaSdHNe1rzS9wbOALaVNEfSmXnc/5ennyvpa4Xf89f5O8yXdHiF73mrpO/0xrZQ2GZ/CNwFjJa0X/5ud+XpNsrjfjXPe76kKYX1tXuOdyZwTGHeO+bfYU7+juPKv4vVqNl3KfrT2A8wlvRwuF1y/zTgw7n7VvIdrMBw4LHcPYl0l+7GwAjSU16PzsPOIT2csTT9Bbn7rcD83P2NwjKGkO7a3zDPt4MKd3IDu5PuqN0Q2Ai4F9g1D3uMfAdu2TTHAxfn7teRHvMxGPg88LVcvgXwl3riyuus9F02AAbn7nFAe+4eD7xESmyDSO+Q+EBZfK8GHgHelPs3IT2ksbN5LiM9z+hVwEzSEwYGk+6YHke6A3kaq+7knwR8P3dPB47M3ccAL+TudYBNCr/zgjyfld8zD9sPmJKHvQr4Vf5d31/6nfN4m1b4LXpzWxhLekLzXoWYbwM2zP0nAF/N3cMK0/0EeE/ungv8Z+4+sxDP90jPuyr9Pus3ex/trx+feQwMj0bEnNw9m7RzduWWiHg+IpaSDmi/zOXzyqafCivfWbCJ0lvb9gNOlDSHdFAZDIzJ48+IiErvNtgXuDYi/hERL5AeDviWLmLcl3TAICIeAB4HticdXA/L43yQ9GgPuhnXusAFkubl+RRfEHZHpPe/rMjrYd+yaV8LLIqIO3OMz0V6xHpX8+yIiFdIj9YYS0qMj0bEQ5GOej+tsj72yXFQWi+ZgG9ImgvcRHrE9+YVpt8vf+4m/cf/OlLCmge8Q9K3JL0lIpZVWX5vbQsAj0d6nwWkl6DtAPwpz2cisFUe9rZ8pjmP9FDOHSVtCgyJiN9XWBczgS9LOgHYKiL+WWX51oV+Vb9r3fZyoXsFsH7uXs6qqsvBnUzzSqH/FVbfbsqfbxOkg9X7I+LB4gBJe5Ieh15JpUdRd6XiNBHxpKSnJe0MHA58sjB+vXF9HlhMevvgq0hnGysXVb7oCvFVev5PZ/Ms/61K67rW5whVGu9DpDPI3SPi30rVk+W/dyneb0bEj9YYIO1OeqbaNyXdGBGn1bDs7m4LlA0TKdEcUTaPwcAPSWfPT0g6NX+vauudiLhC0u2kl5fdIOnjEfG7TuKwKnzmMbA9RqouglVPb63X4QCS9iU9AXQZcAPwmUL98641zOc24BClJ8tuCBwK/KGGaT6Ul7E96T/a0kHqStILqjaNiHm5rDtxbUo6e3gF+Aipiqpkj9zO8CrSevhj2bQPAFtKelNe3sZKjcOdzbOSB4CtJW2b+4+oMt6fSE82hrxeCt9hSU4cb2PVf+3Pk6omS24APlpoTxgpaTNJWwIvRsRPSS+tqvZe7d7aFsrNAvaRtF2exwb59y4lwKdyzB8AiPSWyGU5jtXWhaRtgEci4lxSNd/O3YjHcPIY6M4CPiXpz6R65e54Nk9/Pqvesf11UtXMXKVLJL9ebeKSSK/avYT0FNnbgR9HxN1dTPZDYFCusrgKmBQRpf/cf0Y6kE4rjF93XHkZEyXNIlWJFf8jnklqdJ4PPApcW/ad/kU6oH5P0j2kdpHBXcxzDRHxEul91r9WajB/vMqox5HecX8nKWGUXA60SWonHUgfyPN9mlQVNF/SmRFxI3AFMDOv05+RkssbgDtyldHJwP9WWX6vbAsVvv9SUhvJ1Fz1Ngt4XU4SF5Cq1X5Beq1CyVHAD3KDebFq6nBgfv4urwMuqzceS/xUXbNukDQe+GJEHNTsWFqBpFtJ66O92bFY3/CZh5mZ1c1nHmZmVjefeZiZWd2cPMzMrG5OHmZmVjcnDzMzq5uTh5mZ1e3/A2iC3T+SjGGoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(nread_ovlp, \n",
    "         bins=150, range=(0,150), alpha=0.75);\n",
    "plt.xlabel(\"number of overlap candidates per reads\");\n",
    "plt.ylabel(\"number of reads\");\n",
    "plt.title(\"Distribution of number of overlap candidats per read\")\n",
    "print(\"mean count:\", np.mean(nread_ovlp ))\n",
    "print(\"median count:\", np.median(nread_ovlp) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total reads sampled: 50000\n",
      "total number of read candidates: 2506184\n",
      "number of candidates per read: 50.12368\n"
     ]
    }
   ],
   "source": [
    "nread=len(all_rids[:50000])\n",
    "print(\"total reads sampled:\", nread)\n",
    "print(\"total number of overlap candidates:\", np.sum(nread_ovlp))\n",
    "print(\"number of candidates per read:\", np.sum(nread_ovlp)/nread)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
